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ABSTRACT 


The quasi-static, pressure characteristics of a gas-lubricated 
thrust bearing with shrouded, Rayleigh -step pads are determined for a 
time-varying film thickness. The axial response of the thrust bearing , 
to an axial forcing function or an axial rotor disturbance is investigated 
by treating the gas-film as a spring having nonlinear restoring and 
damping forces. These forces are related to the film thickness by a 

power relation. The nonlinear equation of motion in axial mode is solved 

! 

by the Ritz-Galerkin method as well as the direct, numerical integration. 
Results of the nonlinear response by both methods are compared with the 
response based on the linearized equation. 

Further, the gas-film instability of an infinitely-wide Rayleigh 
step thrust pad is „ determined by solving the transcient Reynolds equation 
coupled with the equation of the motion of the pad. Results show that 
the Rayleigh step geometry is very stable for bearing number, A, up to 50. 
The stability threshold is shown to exist only for ultra high values of 
A 2: 100 where the stability can be achieved by making the mass heavier 


than the critical mass 
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I. INTRODUCTION 


1.1 Introduction 

Among many geometrical profiles which generate hydrodynamic pressure 
in fluid-film bearings, the step geometry is one of the most proficient 
methods to achieve this purpose. This discovery was first made by Lord 
Rayleigh in 1918 (Ref. 1 ), and subsequently there have been many 

contributions (Refs. 2 to 8 ) on the prediction of the performance of 

the Rayleigh-step thrust as well as journal bearings. 

Recently, there has been another significant application of the 
Rayleigh-step bearing in the field of dynamic sealing in advanced air- 
breathing propulsion systems (Ref. 9 ). In this application, a series 

of Rayleigh-step pads is employed on the high-pressure side of a face 
seal in order to maintain a small steady gap (in the order of 0.001") 
between a one-tooth laybrinth and the high- speech rotor surface (Fig. 1 ). 
The flat step is shrouded in order to minimize the side leakage. These 
pads function strictly as hydrodynamic thrust bearings operate in the 
high-ambient pressure , and provide the necessary stiffness to maintain a 
steady gap. The satisfactory operation of this type of seals depends 
critically on the dynamic performance of these thrust pads in the 
presence of an oscillatory force or a disturbance due to rotor run-out 
or rotor unevenness. 

So far, the investigations on shrouded-, Rayleigh-step bearings have 
been restricted to the prediction of the static performance (Ref. 10 ) 
only. The dynamic characteristics of this type of thrust-pads have not 
been given much attention. It is the main objective of this report to 
study the influence of the nonlinear, gas-film, restoring and damping forces 
upon the response of the pad to a given forcing function or disturbance. 


1 
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Specifically, the work reported herein fulfills the following 
objectives : 

1. To develop a gas-film analysis of purely hydrodynamic, 

Rayleigh-step pad to calculate the quasistatic stiffness and 
damping, which depend not only on the operating conditions 
but also on the vibration of the system. 

2. To extend the analysis to the nonlinear axial response of 
the stationary ring due to any external force excitation or 
any disturbances induced by the rotor misalignment or surface 
waviness. 

3. To develop a stability analysis of the infinitely 

wide, single-step pad to explore whether there exists any 
thresholds of stability for the system. 

1.2 Historical Survey 

In 1918, after Reynolds' theory of thick film lubrication became 
generally accepted. Lord Rayleigh (Ref. 1 ) first applied the theory and 
discovered that an optimum profile for the load capacity of a slider 
bearing is a flat step. The method of the calculus of variations was 
used to optimize the shape of slider bearing to an infinitely long, 
incompressible film. From then on, any hydrodynamic bearings composing 
of two sections of parallel surface film are called Rayleigh-step bearings. 

The optimized geometry of a Rayleigh-step slider and its corresponding 
optimized, non-dimensional, load capacity for an incompressible film can 
be found in most textbooks on lubrication (Ref. 11 ) . For compressible 
gas film, the optimized, one-dimensional, Rayleigh-step bearing was analyzed 
by Wylie and Maday in 1970 (Ref. 2), The load capacity "of an optimized step 
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bearing was found to be slightly lower at low bearing number but much 
higher at higher bearing number than the load capacity of an incompressible, 
optimized, slider bearing. 

The problem of a two-dimensional, Rayleigh- step, thrust bearing has 
not received much attention until 1954, C. F. Kettleborough (Ref, 3) 
solved the pressure profile and calculated the load capacity for the step 
thrust bearing by Relaxation methods. In 1955 (Ref.f7)he applied the 
analogy method contributed by A. Kingsbury (Ref, 12) to investigate the 
pressure profile for an oil- lubricated step bearing by an electrolytic tank. 

In 1959, using air as the lubricant, K. C. Kochi (Ref. 6) showed the 
characteristics of an ' inf initely-wide, Rayleigh- step, thrust pad by the 
use of semi-graphical method. He demonstrated that an analytical solution 
to express the pressure profile explicitly is extremely difficult, because 
the Reynolds equation for a compressible film is nonlinear. 

In 1961, J. S. Ausmann (Ref. 7) made certain approximations to 
linearize the Reynolds equation for a compressible film, and obtained a 
series solution of the pressure and load for a self-acting, stepped, sector, 
thrust bearing by the aid of Neumann polynomial. He obtained the numerical 
solutions for the optimized number of sectors and step depth for the maxi- 
mum load carrying capacity. 

Recently, Cheng, Chow and Wilcock (Ref. 8) presented some results 
for shrouded, Rayleigh-step pad used as a flexible face seal. The pres- 
sure generation and static stability of this type of surface profiles 
using as a flexible seal ring was discussed, and the effectiveness of 
hydrodynamic action was confined to the static stiffness characteristics 
of gas film. The influence of the nonlinear gas-film forces and the 
question of gas- film stability of the hydrodynamic, Rayleigh-step pad' 
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has not been investigated. Therefore, to ensure a safe operation of 
Rayleigh- step, thrust pads, it is necessary to conduct a full-scale nonlinear 
study, which is the major objective of this work. 



II. GAS FILM FORCES 


2.1 Statement of the Problem 

The major concern in this section is the determination of the 
time-dependent, pressure distribution within the gas film between two 
annular surfaces containing a series of Rayleigh-step pockets as 
shown in Fig. 1. This problem is formulated within the framework of 
the conventional theory of lubrication for a compressible lubricant. 

The major assumptions commonly used for gas-lubrication theories are: 

1. The pressure across the film is constant 

2. The flow is laminar 

3. Inertia forces are neglected 

4. The film is isothermal 

5. The flow is Newtonian. 

Under these assumptions, the governing equation for a transient, 
continuous, pressure field becomes the well-known, transient, Reynolds 
equation for a compressible fluid. Various analytical and numerical 
methods for the solution to this equation have been outlined by 
Castelli and Peviric (Ref. 13 ). For the present problem, the abrupt 
geometry introduced by the Rayleigh pocket makes it rather difficult to 
solve the pressure equation by any analytical methods. For this reason, 
a numerical method is employed in solving the discretized pressure 
field based on the solution of a system, non-linear algebraic equations 
by the Newton-Raphson procedure. The numerical integration of the 
discretized pressure field gives rise to the time-dependent gas-film 
forces which are required for investigating the nonlinear response 
of this type gas-film to a prescribed forcing function or disturbance 
function. 

- 5 - 
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2.2 Pressure Equation 

The discretized pressure equation is derived by considering a flow 
balance within an element of the gas film as shown in Fig. 2. The mass 
flow rate into the left boundary, AD, and the bottom boundary, AB, are 
designated as q^ and q^. The mass flow rates out of the boundaries, 

DC and BC are, likewise, designated as q^ and q^. In addition, there are 
mass stored in the volume which is designated as q . 

Based on the assumption in 2.1, the velocity is parabolic across 
the film. Integrating the velocity across the film, one may express 
the mass flow rate in terms of the boundary velocities and the pressure 
gradients. These can be written as 


T (urh oh ^ dpi / r i+l r i-l 

i = 1 ? — - 12^ del , . V — 2 

r U h , 3 ^ ,8., - - 0 . , , 

= r„ -L . iL lll\ r 

L P 2 12^ dr J. , . V 2 } i-% 

J 


( 2 . 1 ) 

( 2 . 2 ) 


q 3 


= 'f 0 Jfi* . P h SeI 

! P r> -»« ao r 


L- 2 - izp,T a 0Jj+%(i c 


/ r i+l * r i-l 


(2.3) 


q/. = 


/ r 9 jt . i .. r .. V .i\ r 

L P 2 12n drJ i+%J V 2 / i+% 


(2.4) 


a - iw fo+.l ~ r i-l N, r A+l " e i-l i 

q 5 dt \ 2 J i \ 


l 


(2.5) 


The flow balance requires 


q l + q 2 = q 3 + q 4 + q 5 


( 2.6 > 
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Introducing the following nondimens ional variables: 

e 

P = ^“, H = — , R = ~ , r = u>t , Q = PH , and A = ■— 

*a o p 6 

a 


and using the equation of state p = , Equation (2.6) becomes 

RT 




R... - R. 


A\RH -2- - — . &\ 1 0*1 _ 


L K * 2R 39 A 4 " V™ K ' 2R 90A l V 
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L \3t J. . \ 2 
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/ i \ 2 


Further simplification leads to. 


r /. 
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r * * 1 
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i p + p 

<E 1« - E i-1> L ^ »4*.j 


- “ Of). . ** 

l,j 

Introducing the following finite difference approximations : 

(§A j4 - <q j - • vi>i 


(2.7) 




“ (Q, - Q, ,),/( R, - R, 


^3R/. i . ' x i ~i-l'j 


( 2 . 8 ) 


^i,j“% 


Q. . + Q. . i 
_ 1,3 1,3-1 


2 
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one obtains a system of first order differential equations in time for 
the discretized pressures. 

To obtain the exact, time-dependent, gas- film forces, it is necessary 

to solve this equation simultaneously with the equations of motion of 

the supported mass by an explicit or implicit, numerical procedure for the 

initial -valued problems. This procedure necessitates the calculation of 

the pressure field at each time interval during the transient, and it is 

extremely uneconomical and cumbersome. 

In the case of a high, ambient pressure, the effect produced by the 

term containing ~ becomes insignificant comparing to that by , one may 
OT oT 

3? 

neglect the effect, and thus decoupled the gas-film force calculation 

from the dynamics of the supporting mass. Ignoring the contribution by 
3P 

— 9 the pressure distribution can be solved independently as a function 
OT 

* 3P 

of H and — or H. The assumption of a negligible — effect is made in the 
subsequent analysis. Expanding the terms in Equation (2.7) by using 
relations (2.8) one obtains 


a~Q. . , + a, VQ- • i i • “ (a, + a. + a_ + a,)Q. . + a,Q... . 

3 i, 3-1 1 ^ x i,j-l 5 l-l, j 3 4 5 6 ^ 1,3 6 1+1,3 


+ LH. . , A ‘ a 6j ^i,j + a 4 Q i,j+l " a 2 ^i,j+l 

i, j-i 


= 0 


(2.9) 


where A = H. "t - H. . = i „ 


at the end of the pocket 
otherwise 
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and 
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(H ^» 


(a,) 


(H 3 ), . 
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(H 3 ), 
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x x-1 




Vi,j 4(R i+1 - R i _ 1 )(R i - R i . 1 ) 'i-%,1 


( 2 . 10 ) 


(a c ), , = 


(R i+l * R i ) 


(H ) jjjL , 


Vi,j 4<R i+1 - R. 4 )(R. +1 - R.) 'i+h,i 


(a Q ) . . = AR. H. = AR. ” 
8 1,3 ii i 


2.3 Method of Solution 

The set of nonlinear, algebraic equations to be solved for the pressure 
field, Q. . , is 


I . . = a.Q. . . + a., JQ. . , + 

i,l 3 x i,j-l 1 ^1,1-1 


r a i 


+ a 5 Q i-l,j " (a 3 + a 4 + a 5 + a 6 )Q i,j + Lh7~ 7 A ' + 

i, j“i 


+ a 6 Q i+l,j + a 4 Q i,j+l " a 2^i,j+l * 0 


(2.9) 


where §. . is the implicit function for the node point (i,j). Newton- 

1,1 

Raphson method has been employed here to solve these equations. 
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Based on this method, Q. can be improved by calculating the dif- 

1 » 3 

ference A Q? . which is found from first order approximation of Taylor's 
expansion for §. . . Using Taylor's expansion, eq. (2.9) becomes 



a§. . 

*k.i 




3Q,- 


i, j +1 


A Q 


n 

i, j+1 



(Z.lt) 


From Equations (2. 9) and ( 2.11), one obtains 


- ilj * + ^ AQ ‘‘r + ' 


- ( a* + a 5 + a* ) + I . 


4 


a, 


^ - a 


8 


l 


'V 




Yl- - \ 

+ * (<3« - [ 7 q ^, j 




Cz.n ) 


Since ^ are known from the previous iteration, the set of linear 

simultaneous equations (2,12) are solved for the difference A Q. . , by 

1 j J 

Gaussian Elimination Method. 

The new iterated Q becomes. 


Q 1 ?* 1 = Q n . + A Q. . 

^>1 i»J 


and the procedure is repeated until all A . becomes less than a 
prescribed convergence error. 
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A computer-program RSS-FILM has been written to calculate Q.. for 
a given nondimensionaless , gas- film thickness H and a nondimensionaless, 
velocity H. The nondimensional, pressure distribution can be obtained by 
dividing Q by the nondimensional gas film thickness. H. The horsepower 
required to overcome the friction resistance is calculates as 


T = 


I r 



(2.13) 


and 


HP = 


n x T x N x 60 
63000 


(2.14) 


Also, the -load per pad is obtained by the following integration 


W — p | (P - l)r d0dr 
a J A 


R e_+e T 

2 t* ° » G L 

p r I RdR 
a ° J R. -o 

i 


J 


(p - i)de 


(2.15) 


Figures 3, to 7 show the effects of change in H and H on the pressure 
distributions in a shrouded, Rayleigh-step gas pad. 
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2.4 Approximation for Gas-Film Forces 

The gas-film forces are shown to be dependent upon H, — , and the 
following operating and geometrical parameters 

6p,ujr o 2 

A _ — = Bearing Parameter 


0 G + 


= Step Length Parameter 


2B 

R - R, 
o i 

2uR 

o 

N(R - R.) 
o 1 


= Shroud Width Parameter 


= Length to Width Ratio 


Table 1 lists the gas-film forces for an annular bearing surface 
containing 20 Shrouded-Rayleigh-Step pads, with geometrical dimensions 
shown in Fig. 1 . By plotting data on a log- log scale in Pig. 8 &‘Fig.9 

it is found that they can be fitted with the following function for a wide 
range of dimensionless thickness, H. 


F 





(2.16) 


Tables 2 list the approximate, gas- film force based on Eq. (2.16) and 
the actual values interpolated from Table 1. The errors introduced by 
using the fitted Equation (2.16) are also listed in these tables. 

The values for C^, C 2> n^, and n 2 for the gas-film forces listed in 



Table 1 are found to be 


C x = 2.32 lb 

C 2 = 0.76 lb /(in/sec) 

n 1 = 2.5 

n 2 = 2.5 


The above constants, of course, only apply to the operating and 
geometrical parameters shown in Fig.l , For other parameters, a dif- 
ferent set of C^, C 2 , n^ and n 2 will be required to approximate the 
gas-film forces. 

Letting H q be the equilibrium nondimens ional gas film thickness, 
and defining x as the nondimensionalized displacement from H q , positive for 
a decreasing of H, the gas film force can be alternatively expressed 

as, 


F = 




n„ 


(H - X) 


e dX 

0:6 dT 


(H -X ) 


(2.17) 


where H = H 


X 



III. NONLINEAR AXIAL RESPONSE 


3.1 Mathematical Modeling 

The main problem in this chapter is to determine the dynamic response 
of the stationary ring to any external force excitation or to any disturbances 
produced by the rotor misalignment and by the rotor surface waviness. 

Knowing the detailed motion of the stationary ring with respect to the 
rotor motion, one may calculate time- dependent gap distribution between the 
surfaces. In general, the dynamic response of the stationary ring is 
measured by the axial translation and by the rotations about two mutually 
perpendicular diameters of the ring. However, for a stationary ring with 
a narrow width and a large diameter, the response in the axial mode is 
very weakly coupled with the oscillation in the angular mode. Furthermore, 
the equation governing the angular oscillation is very nearly the same 
as that governing the axial motion. Thus, it is only necessary to 
concentrate the analysis on the non-linear characteristics of the motion 
in the axial mode. This reduces the -problem from a complicated, three- 
degrees-of -freedom dynamical system to a single-degree-of-freedom problem 
for which a more thorough analysis can be afforded. The weak coupling 
between the axial and angular oscillation is demonstrated analytically in 
Reference 1{*. 

Figure 10 shows the mathematical modelling of the non-linear 

vibration of the stationary ring in the axial direction. It consists of 

a stationary ring of mass m subjected to a steady load W D » The back face 

of the ring is flexibly mounted to the frame through a soft spring of 

stiffness K and the front face is supported on a very stiff, non-linear 
s 

gas-film whose restoring force is represented by the power relations 
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formulated in Chapter 2, Eq. (2.17). The major problem here is to 
investigate the motion of the stationary ring for one of the two following 
conditions : 

a, a force excitation, q cos oo^t, acting on the stationary ring. 

b, a prescribed rotor disturbance characterized by a Fourier series 
N 

) e cos nuo,t . 

L, n f 

n=l 


3.2 Equation of Motion 
3.2.1 Force-Excited Motion 

The equation of motion due to a force excitation q cos u>t is 
considered first. Recalling from Eq. (2.16), the force balance on the 
stationary ring gives the following equation. 


A 2 c i 

d x , 1 

m — r + — — r - 

dt 2 h"- 5 


H 


2 dh „ , 

2.5 dt “ V 


k g x — 


q costtJft 


(3.1) 


where k s x is small comparing with W Q and other terms. 


Introducing the following nondimens ional variables 


X =: X / 6 
T - m f t 



m6(<ju*) 
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moo 



m6(oi*) 


where u> is the natural frequency of -the system based on the linearized 
equation of Eq. (3.1). Equation (3.1) becomes 

w 2 X + CX + B [ TT - “111 

(H - X) ZO L (H - X) ZO H Z,:>J 
o o o 

= Q cos T (3.2) 


where 


H « 


H 

o 


- X 


W 

o 


C 


1 



d 

dT 


3.2.2 Displacement-Excited Motion 


N 


cos nu),,t resulted from 
f 


Consider now the disturbance function ) e 

fttl n 

the rotor misalignment, commonly known as the run out, and the non- 
flatness of the rotor- surface. The film thickness H will be perturbed 
by this disturbance, and becomes 


N 

H = H - X + Y E cos n T 
o L, n 

n=l 

€ , 

where E = n ' 5 . 
n 


(3.3) 
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The equation of motion in the absence of the excitation q cos cu^t 


becomes 


m 


A 1 C ! 

dx 1 


C 1 C 2 dh 


dt 2 H 2,5 H 2,5 H 2 * 5 dt 

o 


= 0 


(3.4) 


N 


Substituting Sq. (3.3) into (3.4) anci letting X' = X - Y cos n T , 


one obtains 


n=l 


— 2 ** • 1 

ai X' + CX> 


+ B 


r i 


i i 


2 5 I 2 5 2 

(H q - X') L (H q - X') H q 


.5j 


N 

— / O' cos n T 
L. n 

n=l 


(3.5) 


where 

Q' n = n 2 (3.6) 

Equation (3.5) is identical to (3.2) with the exception that X, and 

N 

Q cos tut are replaced by X' and ^ Q'^ cos nT respectively. Thus, the 

n=l 

solution for the force-excited oscillations is also applicable to the 

displacement-excited motion provided the proper substitutions are made 

for X' and Q' . 

n 

3.3 Linearized Solution 

If the motion of the stationary ring is such that the resulting 
gap variation, H - H^, is only a small fraction of the equilibrium film 
thickness the response can be estimated from the solution of the 
linearized equation about the equilibrium film thickness H q . Linearizing 
Eq. (3.2) for the force-excited motion, one obtains 
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u) 2 X + — X + X = Q cos T (3.7) 

H 

o 

Similarily, linearization of Eq. (3.4) for the displacement -induced 
motion leads to 

-2X4- -~rr X + X = e cos I (3.8) 

0 ) ^ 2.0 

o 

Eqs . (3.7) and (3.8) are clearly the standard, damped vibration equation 
for a single mass, and its solution can be readily written as 


A •= 


2 1/2 

r /n -2,2 C l' 

(1 - (i) ) + " c ! 

H 5 -* 
o 


a = tan 


/ -I 2,„ 

! (1 - id )H o 


(3.9) 


(3.10) 


where X - A cos (T - o') 

3.4 Non-Linear Solution 

Two methods have been employed to obtain the non-linear response 
characterized by the solution to Eqs. (3.2). The first is the' method of 
Galerkin. (Ref. 18) which gives an approximate solution to the nonlinear 
equation. The degree of approximation is governed by the number of terms 
considered in the assumed function in the Galerkm procedure. The second 
method is the direct, step by step, numerical integration using a 
Runge-Kutta procedure. Details of these two methods are given next. 

3.4.1 Method of Galerkin 

The non-linear equation in question, Eq. (3.2), can be represented. 


implicitly as 
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f(X, X, X, T) = 0 


where 

f = w 2 X + G * + 5 - — - Q cos T = 0 (3.11) 

(H - X) H 

o o 

According to the method of Galerkin, one may assume that the unknown 
response X(T) can be represented approximately by a truncated Fourier 
series, 

N 

V 

X = / a cos nT + b sin nT (3.12) 

L> n n 

n=0 

The substitution of (3.12) into the differential equation gives 
arise to the residue function, 

R(T) = f(X, X, X, T) (3.13) 

The R(T) wiLl not vanish unless X(T) exactly satisfy the differential 

equation. 'The Galerkin method provides a set of equations by which one 

can solve for the constants a and b for which the residue function will 

n n 

be made extremely small. These conditions are obtained by requiring that 
the integration of the residue function as weighted by each individual 
Fourier component (cos nT or sin nT) be made equal to zero. Stating 
mathematically, one obtains, 
c 2tt 

J R(T) cos nT dT = 0 

(3.14) 

.2tt 

j R(T) sin nT dT = 0 
o 

for n = 0, 1,2... .N 
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where 


N 


Z — 2 2 

- (ju n (a cos nT + b sin nX) 
n n 


n=0 


N 


V 

B + C i n(-a sin nT +b cos nT) 

n 


> in -a 

U n 
n=0 


N 


H - (a cos nT + b sin nT) - 
i_ o •_ n n J 


f 


n h 


5/2 


n=0 


- Q cos T 


(3.15) 


For N > 1, integration of Equation (3.14) involves definite integrals 
5 

of - ■ th power of the Fourier series. A gallant attempt was made in 
reducing these integrals in some manageable form, but was unsuccesful. 

Thus, the inclusion of any terms beyond N = 1 was not made. 

For N = 1, 

X = a Q + a^ cos T + b ^ sin T (3.16) 

Further simplification is made by representing X with 


X = A cos (T - a;) - A 


where 


A = - a 
o o 


2 21/2 
( a l +b^ ) 


0 ? = 


t -ih 

tan — 
a l 


(3.17) 
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Substituting (3,17) into (3.14), one obtains 


„2tt 


B - CA sin T 


—2 

! - W A cos T H 

° [h_ + A. 1 [ 1 - 


2,5- , -.2. -5 

£ cos t] 

H + A J 


Oj L. 


- Q(cos a cos T - sin o' sin T) J sin T 

J 'cos T 


V dT « 0 


Integrating Eqs. (3.18) and noting the following relations: 


,2tt 


cos T 
sin T 

sin T cos T 


Y dT = 0 


,.2tt 


2 

cos T 
- 

sin T 


dT = rr 


,2tt sm T cos T 


V 1 - rinr cos T ) 

o o 


dT = 0 


,2tt 


sin T 


, x _ & 

( H + A 


cos T ) 

/ 


5/2 


dT = 0 


(3.18) 


(3.19) 


(3.20) 


(3.21) 


(3.22) 


Equation (3.18) reduces to the following algebraic equations 


A (o' 


2-1 Ii - Q cos 0=0 


TT 


(3.23) 


CA 

rr 


I 4+ Q 


sin a 


0 


(3.24) 


l 5 - 2n - 0 


(3.25) 
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where 


— S n cos T dT 
U 1 J o l 


2 ? 2 
— | S 2 sin T dT 

1 ' G 


2H 2 - 5 I 

I 5 = J S 2 dT 

3 o 


r - \ 

2.5 1 + } 


- -Z — -T - COS T 

\ ‘H + A , 

o o 


/n . A _/\ 

, X + — — — cos T | 
\ H + A / 

o o 


D, = (H + A )■ 
1 o o 


Eliminating a from (3.23) and (3.24), one obtains a system of two non- 
linear equations to be solved for A and A q . These two equations are 


2 ,CAI, ,2 


-*-*) + &)' » 2 =° 


G(A,A ) = I_ - 2tt = 0 
o 5 


Using Newton-Raphson procedure, the successive corrections A A, A A q can 

* SF 3F aG 

be expressed in terms of F, G, — , ~rr~, zrr, and —— , evaluated at the 

dA dA dA 3A 

o o 

last interates of A and A . 



o 


3F 

-F 

dA 


o 

-G 

dG 


dA 


(3.33) 



23 - 


A A = 


1 

A 


3F 

BA 

3G 

aA 


-F 

-G 


where A = 


aF 

aA 

ac 

aA 


aF 

SA o 

aG 

3A„ 


TT 


as- = 2/^1 

BA V rt 


A B 

A cu 


ai. 


\ b 

) tt ba 


2 2 

+ 2-^1 


31, 


2 4 ba 

TT O 


oG _ 5I 5 

ba ba 


3G _ ai 5 

aA 3A 
o o 


(3.34) 

(3.35) 

(3.36) 

(3.37) 

(3.38) 

(3.39) 


The differentials of nonlinear integrals can be derived from equations 
(3.36) - (3.28). They are 




■t { 5 jcos 2 T S 3 dl} 


TT 

2 

r 


K i 5 I cos T V T } 


2 "o 


(3.40) 
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Equation (3.40) cont'd. 

■ k ( 5 1 


tj 

-2 


~ = {5 [ sin 2 T cos T S.dT 

SA D„ l J 4 

2 . o 


} 


Tf 


SI 


4 - € { 5 f si " 2 1 S 3 dT } 


ok 


O 2 


TT 


di_ h 2,5 , ,2 

5 _ o 

SA “ I) 


{5 j cos T S 3 dT| 


tt 


SI5 

M 


H 2,5 , ,2 


TT i 5 1 . S W 


where D„ = (H + A ) 
2 o o 


3.5 


and 


(1 


- xs {;} 


H + A 
o o 


cos T 


, J . 3 L -t-J / 

) ( 1 + 


3.5 


H + A 
o o 


cos T 


Being provided with values of A, and A q , or can be solved from Equations 


(3.23) and (3.24) as 


a 


CAI, 


tan 


-2 

BI^ - TT A 0) 


(3.41) 
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A computer program has been written to solve for A, A q and a, from which, 

X and X . can be determined by 
max mm 


X 

max 


- A - 


A 

o 


X . 
rain 



(3.42) 


The integrals 1^, 1^, I,_ and their derivatives with respect to A and 
A q are evaluated numerically. Tables have been prepared for various 
values of A7 (H q + A q ) varying in the range, [0,1] at intervals of 0.01. 

3.4.2 Direct Integration 

The step-by-step numerical integration of Eq. (3.2) is achieved by 
splitting it into two first order equations: 


X = Y 


U) l n 


2.5 


B + CY 


(H o ‘ 


X) 


2.5 


+ Q cos 


T 


1 

J 


( Ref. 16 ) 

The popular Runge-Kutta method has been employed in obtaining the 

A 

solution for X and Y with a given set of initial conditions X , Y . The 
response is represented by the trajectories in the phase space plot (X,Y). 
The response is considered to have reached a steady state if the trajectory 
approaches a limit cycle, which could be a single or multiple- looped cycle. 
A library subroutine at the Vogelback computer center has been used for 
this numerical integration. Fortran listings for the Computer program, 
RSGALN, which calculates A, A^, and by the Galerkin method, and the 
computer program, RSRKIT, which calculates the trajectories in the phase 
space, are included in the Appendices and 
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3.5 Results of Non-Linear Response 

The results of the nonlinear response are presented in two parts. 

The first part is obtained from the method of Galerkin one -term approxi- 
mation, and the second part from the Runge-Kutta direct integration. 

Since the gas-film is an unsymmetrical spring, i.e., the relation between 
the displacement and restoring force is not symmetric with respect to 
the equilibrium position, the amplitude of response during an upstroke 
is different from that during a downs troke. During an upstroke, the 
gas film stiffness is softer, and the amplitude is greater than that during 
a downstroke when the gas-film is considerably stiffer. For this reason, the 
response in the upstroke and downstroke are plotted separately against 
the non-dimensional excitation frequency in Figs. 11 to 19 . 

3.5.1 Results by Method of Galerkin 

Referring to the equation of motion, Eq.(3.2',it is seen that the 
parameters affecting the dynamic response are: 

B = stiffness parameter = 

C = damping parameter 

= static-film parameters = 

Q = forcing intensity parameter = 

a) = frequency parameter — 


m&(ou )' 


c 2 »f 


mm 


7 6 


* 2 
m6(o> ) 


U) 


f/ * 

U) 


The general characteristics of the upper and lower amplitudes as 




2 ; 


a function of excitation frequency are illustrated in the sketch below 


! ^ 

I \ I 

i \1 


/i\ 

/i\ 

/ \ 


The nonlinear response based on one- term Galerkin method is shown as 
solid lines except with a small portion of the unstable oscillations , which 
are shown as dotted lines. The dashed lines show the response predicted 
from the lxnaarized solution. In the region AB, the excitation frequency 
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is smaller than the natural frequency, u) , based on the linear theory, 
and the non-linear solution predicts a smaller lower amplitude but a 
greater upper amplitude comparing to the linear response. As increases 
to the region BC, the non-linear theory yields three possible solutions, 
one along the path BC, one along the path B'C, and another one along the 
path B'C'. The solution along BC is in-phase with the forcing function 
and is the most stable mode of response; the solution along B'C is un- 
stable and only exists mathematically; and the solution along B'C' is 
out-of-phase with the excitation and is less stable than the solution along 
BC. For excitation frequency beyond the region BC, the characteristics 
of the non-linear response are similar to that of the linear response in 

• k 

the region where > u> . In this region, the pad is insensitive to 
the excitation and would not track any disturbance introduced by rotor 
runout or waviness. 

The non-symmetrical nonlinear gas-film produces a response charac- 
teristics which resemble more to the response due to a symmetrical, 
soft, nonlinear spring, for which the resonance occurs at a frequency 
considerably lower than the natural frequency based on the linear theory. 
This correlation is really not surprising, since the mean position of the 
oscillation shifts to the region of softer stiffness, and the nonlinear 
oscillations are dominated by the softer part of the gas-film stiffness. 

Fig. 11 shows typical response curves for the following dimensional 
parameters : 

h = 0.0005 inches 

o 


6 = 0.001 


inche s 
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m =0.2 slug or 6.44 lbm 

= 2.32 lb 

C 2 = 0.76 lb /( in / gee) 

q =1,2, and 5 lb 

The corresponding , non-dimensional parameters are listed in Fig. 11 . 

The inward bending of the resonant peak in the region cu < 1 is. clearly 

visible in all three cases. Fig. 12 shows the effect of increasing the 

mass of the pad from 0.2 slug to 1.0 slug. The increase in mass does 

* 2 

not alter the parameters, B, H o , and Q, since (cu ) is inversely pro- 
portional to m. The only parameter affected by changing of m is the non- 
dimensional damping parameter. A five fold increase in mass is equivalent 
to a times reduction in the effective damping factor C^. The more 
peaky response near the resonance is clearly visible in Fig. 12 when the 
mass is increased by five fold. 

Figs. 13 and 14 shows the effect of increasing the equilibrium film 
thickness from 0.5 to 0.75. The natural frequency is reduced sharply by 
the increase in the film thickness, and the level of response is also 
much greater with a thick film than with a thin film for the same forcing 
function. 

To investigate the effect of damping the value of has been doubled 
and halved from the case shown in Fig. 11 . The curves in Fig. 15 show 
that when the damping is doubled, the response near the resonance is 
considerably suppressed. The opposite effect is introduced if the damping 
is halved as shown in Fig. 16 . 
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3.5.2 Results by Direct Integration 

Both the upper and lower amplitudes obtained by using the step-by- 
step, Runge-Kutta, direct integration are plotted against the excitation 
frequencies in Figs. 17 and 18 . A case of heavy mass, small equilibrium 
film thickness, and large excitation force has been selected to illustrate 
the nonlinear effects. The linear response curve and the approximate 
nonlinear response by Galerkin method are also plotted as dashed and 
dotted lines for comparison. It is seen that the agreement between the 
Runge-Kutta results and the Galerkin results is good near m = 1. This 
clearly shows that even with one term approximation the Galerkin method 
fields a reasonably accurate prediction for the synchronous response. 

For a) < 1, the Runge-Kutta results show a series of superharmonic reso- 
nances at u> approximately equal to 1/2, 1/3, 1/4, etc. The magnitude of 
these superharmonic amplitudes is, of course, governed by the damping 
factor. 

Fig. 19 shows the trajectory in the phase-space plot for condition 
near the second superharmonic resonance. The final limit cycle forms a 
two-loop orbit showing typical characteristics of a superharmonic response. 
Other trajectories at the third and fourth superharmonic resonances are 
also shown in Figs. 20 to 22 . A subharmonic resonance is also found 

for o) ~ 2.0, but the amplitude is small and harmless. The Characteristics of 
the phase space trajectories near W = 1 are plotted in Figs. 23 to 25. 



IV. STABILITY OF AN INF INITELY - WIDE RAYLEIGH- STEP PAD 


4.1 Statement of the Problem 

It is well known in hydrodynamic lubrication that a dynamic system 
involving any fluid-film supports may, under certain conditions, encounter 
detrimental self-excited oscillations commonly known as dynamical instability 
of fluid-film bearings. The gas-film bearings are particularly susceptible 
to this type of instability. The fractional frequency whirl of a 
shaft supported on gas-bearings and the pneumatic hammer in externally 
pressurized gas-bearings are two of the prominent examples of the fluid- 
film instability. For journal bearings, the gas-film instability usually 
occurs xf either the rotating frequency or the supported mass becomes 
large. There have been considerable data available to predict the threshold 
speed or critical mass of the journal bearing; - ' __ .However, 
for gas -lubricated thrust pads, the problem of instability is relatively 
unexplored. Since present trends in gas -bearing are always toward higher 
and higher speeds, it is important to determine whether there exists any 
stability threshold associated with a gas -lubricated, thrust pad. 

This chapter is devoted to the stability analysis of a thrust pad 
with a Rayleigh-step. The geometry of such a thrust pad is shown in 
Fig. 26 . In order not to impose excessive burdens on the analysis, the 
assumption of an infinitely wide pad has been made. Moreover, the motion 
of the pad is assumed to be restricted in the transverse direction only. 

With these two assumptions, the problem is reduced to the stability of a 
single-degree-of-freedom dynamical system with restoring pressures 
governed by a partial differential equation in space and time. 
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4.2 Governing Equations 


The one -dimensional, time-dependent pressure distribution is governed 
by the Reynolds equation, 


f- (ph 3 |e) , + 12(1 

Sx v 3x/ r dx dt 


(4.1) 


where h = 6 + h + e(t) for 0 < x < B, 
o 1 

h = h + e(t) for B, < x < B 

o 1 

e(t) is the upward motion of the pad. 


The boundary conditions for Eq. (4.1) are 

p = p at x = 0 and x = B (4.2) 

3L 

The equation governing the pad motion is, 
d 2 >B 

m — | = j (p - p )dx - W (4.3) 

dt 2 J o a ° 

where m is the mass per unit length of the pad, B the length of the pad, 
and W q the static load imposed on the pad. The dynamic system represented 
by the coupled equations, (4.1) and (4.3) are to be investigated for the 
stability for an equilibrium position. 

4.3 Method of Solution 

The time-dependent, nonlinear, pressure equation, Eq. (4.1), is 
difficult to solve analytically. A numerical approach has been used here 
in solving a set of discretized, time-dependent, pressures along the 
coordinate X. 

The discretization of pressure is achieved by considering the flow 
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balance in an elemental volume within the gas film as shown. in Fig. 26 for 
the i-th element. The flow rates in and out the element are respectively 
and q^, and 


q l = 




1 

2 


q 2 


(. J*L. St + "l + U 2 

\ 12|j, 3x 2 



where = u and = 0 


Considering the flow balance, 

q^ ~ q 2 = rate of mass stored within the elemental volume 

It follows, 

_ (J&L + (jgL. 

\12n 3 J. Vl2p 3# ' , 

lrh% 



(4.4) 


(4.5) 


Using the isothermal relation. 


-^ = R T = constant (4.6) 

p e 

Introducing the non-dimensional parameters. 



X = X /B 


(4.7) 
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Equation (4.7) cont'd 


TT _ Jl _ I* 

6 h 2 - 


6 uBU 


12U.VB 

CT e 2 

Pj 


T « vt 


V « excitation frequency 


and using the following finite approximations for p h , 

\~k i-% 


(§1 


P. - P. - P. - P. , 

i i-I i i-l 


•• ^ ~ X, . 4 X, . 

1 % i i-l i-l 


i (p i + p i-i> 


2 < H i + h-l> 


the discretized pressure equation becomes 


4Tx: ( p iJ - p i 2 ) H. fl - 4^ (h 2 - h-l) ".,1. 


1 /„ 2 


- A (P i+l + P l>< H i + l + V + A (P i + + H i-1> 

= a U L p i <H i + H i-i ) ] 4 x i-i + ai[ p i<Vi + V] AX J ' < 4 - 9 > 


H . - (\l _ , , + H. i’ 

S,1 \ 1+1 1/ 


where 


(4.10) 
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where A. = a .P. , , + d . 

1 0,1 1 + 1,0 0,1 


B. = b .P. + e . 

1 0,1 1,0 0,1 


C. = c ,P. + f . 

1 0,1 1 - 1,0 0,1 


and 


H 


a . = 
o u 1 


s»i_ 


4A X. 

i 


, d . = - A(H _ + H. ) 

0,1 1 + 1,0 1,0 


H 


'o,i 4A X 


s,i-l 


i-1 


, f . = A(H. + H. . ) 

0,1 ' 1,0 1 - 1,0 


b . — - a . - c . e . = d . + f 

0,1 0,1 0,1 0,1 0,1 0,1 


(4.17) 


P, = P =1 
l,o n,o 


(4.18) 


Since the algebraic equations (4.16) are nonlinear, Taylor’s expansion 


is used to reach the following simultaneous equations for A P. , 


i,o' 


A «. = - *.(P. ) 

1 1 1,0 


where A <5. = (A. + a .P.,, )A P... 

1 1 0,1 1 + 1,0 1 + 1,0 


+ (B. + b P. ) A P, 
l o,i i,o i,o 


+ (C. + C .P. ) A P. . 

1 0,1 1,0 1 - 1,0 


(4.19) 


Equations (4 , 19) are inverted directly for successive A P. 

1,0 


convergence is reached. 


until the 
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In equation (4.9), both the discretized pressure P^ and are dependent 
on time. For transient studies, they have to be solved simultaneously with 
the dynamic equation of motion, Eq, (4.3). However, for small oscillations 
and stability analysis, the variation of H and P with time can be con- 
sidered as small perturbations about the equilibrium solution, and P . 
These small perturbed quantities can be expressed as 


H (T,X) = H. (X) + e(T) (4.11) 

1 X 9 O 

P (T,X) = P (X) - e(T) P (X) (4.12) 

x 1,0 x, c 


0 

where e = /6 


!ht| « 1 

The minus sign in (4,12) is due to the fact that an increase in 
thickness leads to a decrease of pressure of the gas film. Substituting 
(4.11) and (4.12) into Equation (4.9), the equilibrium equation and the first 
order perturbed equations can be obtained . Th ey are 


— - — ( p ^ _ p ^ 'i R - — /p ^ _ p ^ \ h 

4A X Vi+l,o i,o/ s,i 4A X x V i,o i-l,o/ s,i-l 


- A (P , + P. ) (H _ + H. ) + A (P. + P. . )(H. + H ) 

1 + 1,0 1,0 1 + 1,0 1,0 1,0 1 - 1,0 1,0 1 - 1,0 


= 0 


(4.13) 


and 


(P_. a . + a .)P.., + (P. b + b , )F , 

i+l,o s,i n,i l+l, c 1,0 s,i n,i i,c 


+ (P. , G . + C .) - d. 
i-l,o s,i n,i i 


3. Be 
e BT 


rp. f. + e.P. 1 
L i,o i i i , c J 


(4.143 



37 - 


where 


a . = 


H . 
S-i. 1 


s, i 8A X. 

x 


H . .. 
s,x-l 


C . = “ 
S " 


, i 8A X. 1 
l-l 


a . = - 7 (H . + H. ) 

n,i 4 1+1,0 i,o 


c = 7 (H. + H. . ) 

n,i 4 i,o i-l,o 


b . = - a . - c , , b . = a . + c 

s,i s,i s,i n,i n,i n,i 


d i * 4TT ( P i+l,c> 2 ' P i,o 2 )( H i + l,o 2 + H i+X,o H i,o + H i,o 2 ) 

- ( P i,„ 2 ' P i-l,o 2 )( H i,o 2 + H i,o^i-l,o + H i-l,o 2 ) 


— (p - P ) 

2 '■ r i+l,o i,o' 


i. = 7 r (H. _ + H. ) A X. + (H + H. - ) A X ~| 

i 4 i_ i+l,o i,o i i,o i-l,o 1-1J 


CT 


f. = - 7 (A X. . + A X.) 
l l l-l l 


(4.15) 


Equations (4.13) are nonlinear simultaneous algebraic equations for 
the pressure distribution P^ They are solved numerically by Newton- 
Raphson method. The boundary conditions are P^( l,n ) = 1. - 
Rewrite equation (4.13) as 


$.(P. )= A.P - + B.P. + C.P. . = 0 

1 1,0 1 1 + 1,0 1 1,0 1 1 - 1,0 


(4.16) 


for i = 2 , . . . ,n - 1 . 
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After the equilibrium pressure distribution is found, Equations 
(4.14) are solved for the first order pressure distribution, P. 

. - - ' 1 5 C 

In Equations (4.14), both P and e(T) are complex quantities, and they 

i,c 

are assumed to be 


P i„c = P i,R + J P i,I 


(4.20) 


and 


e - e e 
o 


3 T 


(4.21) 


The boundary conditions are P. = 0 at entrance and exit, which gives 

i> c 


P, T, = P . = P- = P =0 
1,R n,R 1,1 n,I 


(4.22) 


Substituting Eqs. (4.20) and (4.21) into Eq. (4.14), one obtains a 
system of 2n - 4 simultanous, linear equations, 


A-.- .X + A. .X. + A. - .X. , + A , . X. . . « D. 

x+l,i x+1 r,i i i-l,i i-I i+n-2,i x+n-2 i 


for i= 1, < 


n-2 


(4.23) 


A. .X. + A . X + A X. + A . X = D. 

i-n+ 2,1 i-n+2 i+l,i i+1 i,i i 1 - 1,1 i-l i 


for i = n-1, ...... ,2n-4 

In the above, A and r X - for i = 1, ......... n-2 are defined by, 

j - i+1 


A. - . - P... a . + a . 
i+l,i 3 + 1 , o s, j n,j 


* X i+1 P i+2,R 


A. . — P. b . + b 
1,1 j,o s,j n,j 


, X. = P 

i i+l,R. 


(4.24) 


A. , « P. , C . + C . 

i-l, i 3 - 1,0 s,j n,j 


X. = P 
i-l i,R 


A • , rs . "" 6 
i+n-2,\ j 


’ X i+n-2 P i+1, 1 


D. = d. 
i 3 
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and for i = n - 1, . , . , 2n - 4, they are defined by, 
j «= i-iri-3 


A. = P., , a . + a . 
i+l,i j+l,o s,j n,j 

A. . = P. b . + b . 

i,i 3,0 s,j n,j 

i-l,i 3 -l,o s,j n, j 

A ~ . = - e . 
i-n+2,1 3 

D. = f.P. 
i 3 3 >° 


X i+1 ‘ i-n+4,1 
X. = p 

i i-rt+3 ,T 

x — p 
i-l ~ i-n+2,1 

X i-n+2 = P i-n+3,R 


(4.25) 


The boundary conditions become 


X 

o 



n-1 



for equations with i = 1, 
n - 2, respectively 


(4.26) 


and 


X n-2 " *1,1 - ° 


X„ = P = 0 
2n-3 n,I 


for equations with i = n - 1, 
2n - 4 respectively 


(4.27) 


These simultaneous equations (4.23) can now be solved for [X] vector 

by direct matrix inversion. A computer program has been written to solve 

for P. and P. for different real values of v, and the Fortran listing 
1,0 i,c 

of this program is included in Appendix D. 

Once the real and imaginary part of the gas film pressure are determined, 

the integration of P. and P. gives respectively the in-phase and 

X y R X, X 
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out-of -phase bearing forces. The in-phase force P dX, can be interpreted 

«J _ R -j 

as the stiffness of the film, whereas the out-of-phase, j P dX, represents 
the damping factor of the film. It should be emphasized that both the 
in-phase and out-of-phase perturbed pressure are dependent upon a, or 
the excitation frequency v of the mass. The frequency-dependent charac- 
teristics of the gas-film reactions is a direct consequence of the inclusion 
3P 

of the term — . The use of these frequency- dependent bearing forces in 
determining the dynamic stability of the gas-film and pad system is des- 
cribed in the next section. 


4.4 Stability Criterion 

The stability of the gas-film and pad is gov'emed by the equation of 
motion, Eq. (4.3), which, in its non-dimensional form, appears as 


-.2 p B 1 W 

— | = “^2 f (P ‘ 1 ^ dX H 

dT m6v o m6v 


(4.28) 


Recalling the pressure is the summation of the equilibrium pressure, P , 
and the dynamic pressure, - e? c , one obtains 

p B ^1 W p B 1 

| (P - 1) dX I P dX (4.29) 

m5v o m6v m6v d o C 


It follows that 

,2 .1 

) — f + e J P c dX = 0 (4.30) 

dT 0 o 

Mathematically, Eq. (4.30) represents a free oscillation problem which 
contains stiffness and damping factors depend ing upon the frequency of 
the oscillations. A direct approach in determining of the stability of this 


/mSv 

\p B 
a 
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single-degree of freedom problem is to look for the eigenvalue of this 

system. If the real part of the eigenvalue is negative, the system is 

stable; otherwise it is unstable. If the eigenvalue is purely imaginary, 

then the system is at its threshold of instability. 

Alternatively, one can also determine the stability thresholds by 

assuming that the eigenvalue is a purely imaginary number and inquire what 

m6 2 

would be the mass parameter, — for a purely imaginary eigenvalue. 

Let the eigenvalue be represented by X + jv, and for a pure imaginary 
eigenvalue, X = 0. It follows that 

e = ® 0 e ^ Vt - e 0 e ^ T - (4.31) 

Substituting Eq. (4.27) into Eq. (4.26), and separating the real from the 
imaginary part, one obtains 

J Pj (v) dX = 0 (4.32) 

0 o 



(4.33) 


where P.j.(v) and P^(v) are solutions of Equations (4.23) for a given value 

of v. The pure imaginary eigenvalue v may be determined by evaluating the 

..1 

integral, I P (v)dX , for various values of v until the integral changes 
J o 

its sign. The exact eigen frequency v may be calculated by a linear 
interpolation. Once the eigen frequency is found, the critical mass at the 
threshold of instability can be determined by Eq. (4.28), and 
P » -I 

*ct " ~2 , V V)dX 

6v o 


(4.34; 
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Equation (4.29) predicts a quantitative value of the critical mass, 
but does not furnish any information on which side the stable region lies. 

To determine the region of stability, one may use the criterion developed 
by Malanowski and Pan (ref,. 17). Their criterion can be stated in the 
following manner. 

Stability Criterion - The system consisting of a thrust bearing of mass, m, 

with Rayleigh Shrouded-Step Seal is in a state of self-sustained oscillation 

at frequency v when and only when the out-of-phase component of the 

° ,1 

bearing reaction, J Pj(v)dX, vanishes, and when the mass, m Q has the 

o 

critical magnitude to be in resonance with the in-phase bearing reaction. 
Equation (4.28). The system will become unstable if the mass exceeds the 
critical value provided the out-of-phase bearing reaction increases with 
the frequency in the algebraic sense in the neighborhood of the critical 
frequency V 0 and conversely. 

4.5 Results 

The steady state pressure distribution, P q , has been calculated 
numerically for the following parameters: 

B l/_ = 0.5, and 0.75 

15 

H ~ = 0.5 and 0.75 

o 

A - « 8.4, 42, and 100 

Pa 6 

The resulting pressure curves are plotted in Fig. 27 > and they are in 

excellent agreement with the analytical solution provided by Kochi ( -Ref. 6) 
This comparison confirms the accuracy of the present numerical solution 
of the steady state pressure distribution,?^ 
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For each steady-state pressure distribution, the dynamic pressure 
distributions, P R and P are calculated for a series value of ct. Figs. 28 
to 35 show a typical series of the dynamic pressure profiles for B^/ B = 0.5 

and 0.75 , H = 0.5, and A = 42. For small excitation frequencies, the 

real part of the dynamic pressure is dominated by the bearing parameter A, 
and the profile is similar to the static pressure distribution, and is 
relatively independent of c. As ct increases, the pressure distribution, 

P R , becomes slightly wavy at both edges. The waviness penetrates deeper as 
ct further increases. As cj approaches infinity, the effect of A disappears, 
and P R approaches the asymptotic solution, 

P 

(p ) - _2 

v R ct - « H q (4.35) 

The imaginary part of the dynamic pressure takes a wavy pattern even for 
small values of ct. For ct approaches infinity, the values for P^ vanishes 
throughout the entire region. 

The in-phase and out-of -phase forces are plotted as a function of 
ct in Fig. 36 , and also listed in Tables 3 to 5. Xt is seen that for» small 

or moderate values of A, the out-of-phase force never becomes negative. 

This indicates that for nearly incompressible cases, there exists no stability 
threshold, and all equilibrium solutions are stable. As A becomes ex- 
tremely large, the out-of-phase force does become negative at a fairly 
high value of ct. This crossing-over of zero line indicates that for a 
highly compressible film, there does exist a stability threshold, and the 
gas film will exhibit a self-excited oscillation at a fairly high frequency. 
Moreover, the criterion in Reference 16 suggests that the stable region 
lies in the area where the mass of the pad is greater than the critical 
mass calculated according to Eq. (4.34). 



V. SUMMARY OF RESULTS 


1. The gas-film restoring forces in a Rayleigh- Shrouded- Step thrust 
pad can be determined numerically by solving the discretized time -dependent 
Reynolds equation with irregular grid-spacings to account for any abrupt 
changes of pressure at the step and at the exit edge. 

2. For conditions of high ambient pressure, for which the term of 

h can be neglected in comparison with other terms at the right side of 
ot 

the Reynolds equation, the gas- film force is found to be approximately 
inversely proportional to nth power of the film thickness and directly 
proportional to the squeeze-film velocity. The exact value of n is a 
function of the step geometry. In general, n lie between 2 and 3. 

3. The axial, non-linear response of the Rayleigh-Shrouded-Step 
pad to a sinusoidal, axial forcing function or a sinusoidal disturbance 
due to the rotor misalignment or surface waviness can be determined by 
one of the following two methods : 

a. By assuming the response to be a truncated Fourier series 
in multiples of the excitation frequency, the Ritz-Galerkrn 
procedure can be employed to predict the non-linear behavior 
of the pad motion. 

b. By integrating directly the equation of motion of the thrust 
pad using a step-by-step, numerical routine, the Runge-Kutta 
procedure. 

4. Results obtained by using the Ritz-Galerkin method with the 
first harmonic terms show considerably departures from the linear response 
curve as the frequency approaches the resonance based on the linear theory. 
The asymmetric spring characteristics of the gas film result into a non- 
linear response similar to that caused by a symmetric, soft, non-linear 
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spring. The resonating peak bends inward and occurs at a frequency less 
than the resonating frequency based on the linear theory. The peak can 
be suppressed by decreasing the mass, increasing the damping, and increase 
the stiffness. 

5. Results obtained by the step-by-step direct integration confirms 
the approximate solution in the vicinity of the resonance. The direct 
integration also predicts a number of additional peaks at frequencies less 
than the resonating frequency known as the superharmonic resonance. 

6. The gas-film instability of an inf initely-wlde Rayleigh step 
thrust pad canbe determined by solving the complete, time dependent, 
Reynolds equation coupled with the equation of the motion of the pad. 
Results show that for bearing numbers. A, up to 50, the Rayleigh step 
geometry is very stable, and no stability threshold has been discovered. 
For ultra high values of A £ 100, a stability threshold is shown to exist, 
and the stability can be achieved by making the mass heavier - than the 


critical mass. 
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TABLE 1 


GAS FILM FORCES (lb f ) 

d = 6.46 ins d. = 5.96 ins 
o 1 

6 = 0.001 in (step) 
uj = 277 rev. /sec. 

P a F s 315 psia.: Ambient Pressure 

a o r 


' — H . 
Nmn 

dh/dK^ 

2.0 

1.5 

uo 

0.75 

0.50 

0.30 

0.20 

-1 in/sec. 

0.50635 

1.10134 

3.2606 

6.9965 

19.949 

65 . 138 

139.729 

-0.5 

0.45104 

0.97913 

2.8878 

6.1765 

17.499 

56.041 

115.531 

-0.25 

0.42338 

0.91801 

2.7015 

5.7667 

16.275 

51.529 

103.734 

+0.25 

0.36806 

0.79580 

2.3290 

4.9477 

13.833 

42.576 

80.750 

+9.5 

0.34040 

0.73470 

2.1428 

4.5384 

12.615 

38.135 

69.559 

+1 

0.28509 

0.61250 

1.7704 

3.7204 

10.182 

29.325 

47.778 

0 

0.3957 

0.8569 

2.515 

5.357 

15.06 

47.124 

92.8 

F - F 
gl g2_ 

0.1106 

0.2444 

0.7450 

1.6380 

4.8836 

17.907 

45.975 


fi 2 - 

lb - sec 
p, * 6 x 10 3 ^ — 

in 

A = 8.30076 
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TABLE 2 


LOAD (lb f ) CALCULATED AND ERROR OCCURRED (%) 


\H . 

.\ min. ... 

h in/secjv 

2.0 

1.5 

1.0 

0.75 

..50 .30 

0.20 

-1 


0.544 

1.116 

3.08 

6.33 

17.403 62.7 

171.2 



+ 1.5% 

+ 1.45% 

- 5.53% 

- 9.5% 

- 12.5% - 3.7% 

+ 22.5% 

-0.5 


0.477 

0.978 

2.70 

5.55 

15.28 54.95 

150.1 



+ 5.78% 

- 0.1% 

- 6.5% 

- 10.1% 

- 6.29% + 6.65% + 30% 

-0.25 


0.4435 

0.909 

2.51 

5.16 

14.205 51.075 

139.05 



+ 4.78% 

- 0.98% 

- 7.07% 

- 10.5% 

- 12.3% - 0.89% 

+ 33% 

+0.25 


0.3765 

0.771 

2.13 

4.38 

12.055 43.325 

118.95 



+ 2.34% 

- 3.01% 

- 8.2% 

- 11.5% 

- 12.9% + 1.76% 

47.4% 

0.5 


0.344 

0.703 

1.94 

3.99 

10.98 39.45 

107.9 



+ 1.18% 

- 4.22% 

- 9.35% 

- 11.8% 

- 12.9% + 3.46% 

55.2% 

1 


0.276 

0.564 

1.56 

3.21 

8.83 31.7 

86.8 



- 3.15% 

- 7.85% 

- 11.85% 

- 13.7% 

- 15% + 8.2% 

+ 82% 
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TABLE 3 


CT 


20 

50 

100 

200 

300 

450 


DYNAMIC BEARING REACTION FOR 
Bj/B = .75 H 2 = 0.5 


A = 8.4 


,1 

r 1 


J V* 

o 

o 

13373 

1.1836 

21003 

1.2694 

27600 

1.4097 

24248 

1.5930 

18121 

1.6590 

12761 

1.6887 


A = 42 


1 

.005975 

10 

.057536 

50 

.11984 

100 

.020667 

200 

.031105 

300 

.26865 

500 

.37916 

2 x 10 3 

.079279 

4 

10 

.045201 


1.5018 

1.5144 

1.6844 

1.7192 

1.5013 

1.5269 

1.8392 

2.0095 

2.0296 


A = 100 


10 

.05897 

50 

.28043 

100 

.47831 

200 

.47977 

300 

.13904 

500 

-.12737 


1.0295 

1.1019 

1.3014 

1.7884 

1.9829 

1.5789 


- 50 



TABLE 4 


<7 


10 

50 

100 

200 

300 

500 


1 

10 

50 

100 

200 

300 

500 


1 

10 

50 

100 

200 

300 

500 


DYNAMIC BEARING REACTION FOR 
Bj/B = .75 H 2 = 0.75 

A » 8.4 

r 1 

\ V* 


I r a“ 

O 


. 15200 
.16382 
.15032 
. 14273 
. 12283 
.08910 


.55179 

.78737 

.84853 

.92439 

.96835 

1.0045 


.009264 

.090506 

.27479 

.17722 

.052296 

.10195 

. 15700 


A = 100 

.004522 

.045139 

.21611 

.37708 

.42895 

.22267 

-.036699 


.61777 

.63251 

.86727 

1.0606 

1.0238 

.99666 

1.0797 


.46999 

.47216 

.52320 

.66548 

1.0325 

1.2315 

1.0725 
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TABLE 5 


a 

1 

10 

50 

100 

200 

300 

500 

800 


DYNAMIC BEARING REACTION FOR 
B^B = 0.5 H 2 = 0.5 


A = 8.4 


,1 

r 1 

J V* 

o 

J o R 

.04398 

.9242 

.40067 

1.0473 

.65404 

1.9523 

.34743 

2.2234 

.17872 

2.2381 

.14173 

2.2543 

.10469 

2.2785 

.07845 

2.2960 


A * 42 


1 

.00519 

10 

.05146 

50 

.22100 

100 

.38003 

200 

.66419 

300 

.65252 

500 

.22427 


2.3313 

2.3349 

2.3947 

2.4791 

2.7731 

3.1907 

3.3674 


1 

.00727 

10 

.07261 

50 

.34848 

100 

.61342 

200 

.73899 

300 

.49508 

500 

.29231 


= 100 

1.9138 

1.9170 

1.9943 

2.2089 

2.7584 

3.0561 

2.8927 
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r 0 s 3.23 in 
r ( =2.98 in 

8=0.001 in (step) 
< j )~ 1740 rad/sec 


08 


P 0 =P 0 =3I5 psia 
A= 8.30076 


e G = 10.9° 
0 L = 5.32° 



Fig.l The geometry of a shrouded, Rayleigh-step 
thrust bearing 




54 - 


if 



Fig-2 Flow balance around a typical grid point 
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Fig. 5 Pressure distribution for h*- Ups. l 
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Fig. 7 CONTOUR. MAP FOR PRESSURE DISTRIBUTIONS 
for h =-l in/sec, H^jfO.5 













63 



OJ 


0.2 

At 


0.3 



Fig.12 ' Nonlinear response for H<r0.5, m=l slug 



(J=976 rad/sec 


Q=0.I27 

S — — - -- 

Q= 0.06 34 
^ Q=0.03I7 
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Fig.l 5 Nonlinear response for 


H 0 = 0. 5, m* 02 slug,C «1.52 It^n^ec) 
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o/=!984 rod/ sec 



Fig. 16 Nonlinear response for H* e 0.5 
m*0. 25 slug C 2 =038 Ibf/i ps 


0.8 



Fig. 17 Comparison of the approximate c 

of re$pons< 




Q s 0.076 

Ho 0.5 
m= I slug 

-I- Runge-Kutta 

— Galerlfin Method 

— Linearized 

Prediction 

* 

cu* 887 rad/sec 



0*0.076 

Ho 05 
m=l slug 

+ Runge Kutta 

— Golerkin Method 

— Linearized 
Prediction 

*887 rad/sec 


Fig.18 Comparison of the approximate and exact downward amplitude 

of response 
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harmonics 
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-0.54 


Fig.2l Phase plot with forth harmonics 
a>=0.248 
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Fig.22 Phase plot with 2nd order 
subharmonics 

a) = 2.096 
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< uu t -o.jssn 
Co~0.1oc 


Fig.2i3 Phase plot with limiting cycle of 
large amplitude 
cj =0.732 
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Fig.24 Phase plot with limiting cycle 
of small amplitude 
5=0.745 




Fig.25 Phase plot at 
aj =|,0 



Fig.26 Flow balance and geometry of a infinitely 
wide Rayleigh step pad 
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Fig.27 Pressure profile for B,/B=0.5 t H 2 =0.5 
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e for J{=42, 













Fig. 36 Variation of dynamic bearing forces with the excitation 
frequency, CT for B,/B = 0-75, A=42, H 2 =0-5 
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PROG RAM T'SEA'CI r NPtTT'VOtJTP U TTT'AP E"5 =TN PUT'? TAPt 6 = 0 CTT P U TT ~* " 

RAYTH 

DlWENSI'0N' DR'(17TVDTHr33 ) fHMTNA ( TO'} jRPSA( IOTjPAAX l'O )TVTS AT10 ) >TH( 33' 
1) » R ( 17 ) » FRC (.17 ) » H ( 1 7 » 33 ) *A1( 17*33) »A2(17»33) * A3 ( 1 7 » 3 3 ) » A4 ( 1 7 » 33 ) s 
' '2 A5T17V3'31 rA5TT7 V S3 ) »Ag( I7»33 ) *6(17) »C ( 1 7 )~» A ( IT* 3 3 )' » F ( XT) »QSMAT17» 
317) »E( 17,17 ,33) »G( 17,33) *DQ( 17,33) ,P (17* 33) * Q ( 17 , 3 3 ) ,QQQ ( 33 ) » 
tpp c r7'r , "ftoutt roy ,tdo tctot — — - ™- - 

1 FORMAT ( 72H 

" 1 ' ” F ~ ~ 

2 FORMAT (1615) 

“3 FORMATTS'F l'O «T) — ' ~ 

4-O-MAT (4-5 »5X*7ET0.4/ (8E10C4DD 

~5 T0'R''M'A'TT'6'X7 AH'l A^TT 6X7'4’FTN PATT* 9X i * 1 FTM f VX TRW* 8 XT’ TFT IT* 3 XTTR J HT7 X ~ 
1 3H I HH» 7X » 3PUHH, /) 

'”6 F CFRM'ATTSTrcr* / 1 — 

7 FORMAT ( AX » 6HLKOUNT»5X 5HNDIAG»6X 4HIRRG* 6X * 4HNPRE* / ) 

8 TORMAT(7F5iH “OUTSIDE ’DIAMETER ( INCHEST" " ~ 

1,E12.5/) 

‘ 9 FORMATTSSH TNSTD'E"DTAMETER( INCHEST — — , 

1 E12 • 5/ ) 

10 FORMATT5'5H THt ANGLE" EXTENDI N6 ' THE 'POCKET' REGTONI DEG - . ) ~ > 

1 E12.5/) 

'TT FORMAT ('5*5 H THE ANGLE" “EXTENDI NG THE LAND 'REG ION (’DEG.) * - ■ ' *' 

1 E12 • 5/ ) 

~T2 F0RMATC5TH STEP" DEPTHTT N CHE'S) ' ' »" 

1 E12.5/) 

T3 FO-MATT5S4-- 'OOTE ; - Tl DTH OF" THE SH-0UDT1NCHE'S ) ------ -- - , 

4-124-5/D 

"T4 -FORMATT55H*TNMER"WTDTFT O'F ' THE~SH-OOD'(TNCRES ) » 

1 E12.5/) 

"T5 TOWATTSTH' 'CONVERGE RTF' 'ERROR' — , 

1 E12.5/J 

"T6 FORMATT5AH GRI D SPACING S IN ” T HE RADIAL DTRECTTON CIWCHES ) ‘ ' " / ) 

17 FORMAT ( 54H GRID SPACINGS IN THE CIRCUMFERENTIAL D I RECT I ON < DEG ) /) 

T8 FORMAT ("5 5 H “MINIMUM - FTLM THTCKNESS (INCHEST " "■ ' , 

1 E 12 • 5 / ) 

T9 "FORM ATTpSH 'REVOUUTTOITS'" PER" SECOND" — ““ — ~ T 

1 E12.5/) 

'TO ' 0 — — TTTTH' T El 2. '57 — “ * ” — 

2 -4-2S -E0-S4-4- , 

TET2'.57T““ — - — 

22 FORMAT { 55H VISCOSITY ( LB-SEC/4-N$$2 ) , 

T - ET2T5T) — — — • 

23 FORMAT ( 55H LAMDA ( 6. *VI S$6. 28SRPSS ( RO/HMI N ) **2/PA ) » 

r ET2C57D' " ' ----- 

24 FORMAT ( 12H RBU)*TH(J) ) 

TT ' fo-m‘ATT"2'5H"rl;h-thi:hthrrthth',hbh- “7“ i •“ — -■ 

2 — 1TU15+ A1U4 -T JDTET-C D 

27' FO- MATT ATI-FiT)' FEOFF D'rCTTDTAT+'TI r=r r ) TAjf T~» TT»AT 4-TTFTT"' “TT0X7THI = 

1 ♦I5*10X»2 HJ=,I 5) * 

-'28 FORM AT(-~2'5TT “MATRTX" TSTT NGULAR AT J- - 1 3 »T6H , CASE ABANDONED". /TH'l ) 

29 FORMAT ( 29H0F INAL PRESSURE DISTRIBUTION. //) 

TO 'FOR MATT" 7/T8H TA'SE~ CONVERGES TO F97ST' 

1 I4H 


AFTER 1 3 ♦ 1 1H ITERATI 



31 FORMAT ( / 10( 1X»F11,7) > 

3’3 FORMATT55H FI'L'M' TR1 CKNESS AT O.D' (INCHES! * ■ " 

1 E 12 • 5 / ) 

" “ 3*4 FORMAT T7'5'5H' TOTATT XO'AD'fLBSl * " 

1,E12.5/ ) 

J5~ FORMAT(T5Tr'DTME71TI ONUES'S' T0~AD=L70AD7TATTEA^P7Y} 

1 E 12 • 5 / ) 

36‘ 'FORMATT5‘5T H(JRSEPOWE-“[:OSS “ ' ' ' — — — i • 

1 E12.5/»1H1) 

3"7" FORMAT.("5'5H ST ATTC "ST3UEF2.E 'PI LM \7ET_0 (TN7~SEC ) ' — , ~ 

1 E12.5/ ) ' 

- ~3T‘FORMAT~n 57S3K1 T8ET0V4T) * ‘ — ~ ~ ' ~~ 

3 +0-MATU55+ D+MEN5+-NLESS STAT+C F+LM '.VELOC+TY W+TH HO = STEP DEPTH t 

1 "E F2T57T ' ” 

N R ~5 

READ! NR * 1 ) 

’ ' 9TT ~RE7CD rWRT21 DKST? NTP'A'D *~MTN »TH» J H » I HH T“J HfT» LXOU1TT , KDTAGTTRRG *T1PRE TNSYM 
1 »ND1 

RE AD! NRTTTDCTTDT ,THG> TBL7ST EPD , WO , W‘I TERROR ' - ‘ 

NN=N-1 

WRITE TNW » 1 J 

I'FTI RRTTTMEVITG'O' TO'" 9T — 

READ ( NR » 3 ) ( DR( I } ♦ 1 = 1 »MM) 

RFATJ { -R 1 R-,“3-)TDTH-(-JT»“0 = r 7NN T ~ 

WRITE ( NW» 16 ) 

WTTETmflT2Trr(W(TT, r^TYMMT — - - — 

WRITE (NW, 17) 

WRTT ETTWr2TrrETH rj)TJ=TTNN'r - - --■ 

TWD = I DO-D I ) *0 4 5 

' ■ ITTMS7M7EQTIT 7 WD'=TWD*tf; 5 ‘ * 

tlg=thg+thl 

D’Cr _ 70"0 T=X »15fM --- • - 

700 DR( I )=DR( I >*TWD 

DO 70 f TTiViHR ~ 

701 DTH(J)=DTH( J)*TLG 

GO ~ T 0" ~9 5 ~ ~~~ “ 

_91 FH^H" 1 

IHH1= 1HH— 1 

; — '■* 

bihi=m-ihh 

IHTTIH-I -~ i ' 

AIH1=IH1 

rrrwsYWTnsTTj ~go~io~jtc — ^ j — 

IF C IH .EQ* 1 ) GO TO 709 

D - 0 - 9:2 - — r^ITITHT — — * 

•92 DR(I)=W'I/AIH1 

~~7XT9 ' U0~97 T='nT»TH HI * ' ‘ — * - 

97 DRU )=( (DO-DI ) *• 5-WI— WO ) /AIHHl 

'GCrTTT^O — — 

710 F I HH1= I HH 1 



rEMP=T {'DO“DT)^0“.'Z5-WOTTF ITCH I~ 

DO 711 I = 1 » I HH1 

711 D'RTI )' = TEHP * ----- - 

720 DO 98 I = I HH »MM 

"98 DR7 I )'=W07BH-T1 ' ~ " 

JHH1=JHH-1 

AJ'fl'Ffl='JRTTl ’ 

B JH=N- JHH 

"D'O WJ = r, JHH1"" ■ ' 

93 DTH( J) =THG/AJHH1 

DD *WJ5 CIHKiW 

99 DTH( J )=THl/BJH 

~WTtwronjE* 

READ ( NR * A ) NV I SM » l V ISA ( I ) * I = 1 * NVISM ), 

RFA'DTN RTATTKPT MITRF'S'AT IT*7 = l‘iNRPSMJ 

READ< NR* A ) NPAK , ( PAA ( I > » I=1»NPAM) 

R EA D CNR' ,“A )"N HMW * 1 HM7 N AT‘I )' ,"1=1 ,'NHMM ) " 
WRITE INPUT DATA 

READ' ( NRT31 rRUOTTTTTr^Tl NRIWT 

READ (NR, A) NHDTM * ( HDOT ( I ) » I = 1 »NHDTM ) 

WRTT ETNWT5 ) — - - - 

WR I TE (NW, 6) LAST, NPAD»M,N» I H»JH, I 

■ ' WRITE (NWT71” 

WR I TE ( NW » 6 ) LK.0UNT »NDI AG » I RRG » NPRE 

WRTT El NW8TDCT ‘ 

WR I TE ( NW » 9 ) D I 

WRTTETTIWiTO TTHG “ ‘ 

WRITE (NW, 11 JTHl 

W RTTETRW El 2 ESTEP D 

WRITE (NW, 13 ) WO 

' WRTTETHWT TATWT “ ' 

WR I TE ( NW , 15 ) ERROR 

* ”1 F*{ N DTA'GOT 3TTT3 1T73 0 — * ~ 

730 WRITE(NW,16) 

— WRTTETm*T2 OTTORTI T»T-TTMM) — 

WR I TE ( NW , 17 } 

W'RTTETR'WT 2 OTCDTHTOTj J= r » NNT ~ ” 

731 CONTINUE 

DO - 1 0 0 0 ‘"N VI 'S= IT NVTS M — ~ 

VIS=VISA(NVIS) 

- - D "CTT 0 O'O' * 'N R'P S = X » NR'P'S M ~ 

RPS=RPSA(NRPS) 

DO" TOO'0"T1FA^T, MPA'M “ 

PA=PAA(NPA) 

D0~9 _ 6 _ "T=1TM" " 

DO 96 J=1 »N 

— 9'6 QlT ,T)‘ =1T ** 

DO 1000 N+DT31TN+DTM 

DO "T00'0'"NHr=T= I'TNRWM ' 

HMIN=HMINA(NHM) 

CONE=THOlXn NHH r/TTHTN ' 

CONE1=CONE-1. 

WRTTF TWWV3T) FTOTTTTNHM") "* “ 

WRITE(NW,18 )HMIN 



“ ¥R"rTE(WVl^TRP'S - 

WR I TE ( NW , 2 1 ) PA 

WRTTE’t NW, 22 fVIS" " ~ 

0-+TE UN0T3 D +DOTUN+DTD 

' " pr=3“.T2as — 

RO=DO/2.0 

'RT=DT7Ti 

PL AM=6 • *V I S#RPS*2 • *P I * ( { RO/HMIN )**2« }/PA 
CHDO T=H DDTfffflUT )7T2 . *P I * R P S*HM'I N > 
CriDT=HDOT (NHDT )/( 2.*PI*RPS*STEPD> 
WRTTETnW'* 39T ”CHDT' ' ~ 

step=stepd/hmin 

WRTTEITTWT 2TTTP1 LAW r--- 

STEPH=STEP*0.5 

C C=XT« 0 T 74- 5‘3'2~9 

GENERATE COORDINATES 

Tflm=OV 

NN=N-1 

DO ”10 0 ■J=TTTTN' 

100 TH( J+1)=TH( JH-DTHl J}*CC 

” R ( 1 ) = D r/'D'O ' - - -- -- - 

MM=M-1 

D'0'TD5"T = r,MM' " 

105 R< I+1)=R{ I )+DR( I >/RO 

DTDO”=DT7D'0 “ ’ ■ 

OD I D0=1 • ~D I DO 

’L6~= 1 ” ’ “ ' •" — 

K.8- 1 

AUHr=T.’ - 

109 DO 107 J= 1 » N 

~TT7 H'CC'6 CJT^A'DNr _ ' ‘ — 

IF (K8.EQ.2) GO TO 108 

L6=M 

K8 = 2 

ADMT=CONE 

GO TO 109 

~ roe ‘ w tost i =”2 rm — 

RS9=R ( I ) 

D0‘1W~7J=1 ;TT ” “ 

106 H( I » J)=l*+CONEl*( ( RS9-DIDO)/ODIDO)**2 
IFTNDTATT) * IT4 a 1T4TTT0 

110 WRITE ( N W » 24 ) 

W RTT ETIWT2~0XrRTT7TT =TTfi~) ~ " ~ 

WRITE(NW»20> (TH( J) ,J=1,M) 

■ GENERATE ATT 1 ,73') 'T'CT ~A“6 CTTS) — 

DO 112 I = 1 ,M 

TT2' IffRTTE TRWT2T) ) “ THT IT JJVJ » 1 *RT) 

114 DO 140 J=2 , NN 

D'ZT^rmO+1 r=TrHT3=T3 ■ 

DZ2=2.*DZ1 

D~Z 3'~ DZT2' 1 * r CTHTJT— THT J— "ITT 

DZ3= 1 . /DZ3 

DT^=A'.V'('D27^TTH ?c;vrr-THl'J ) ) 1 ~~ 

DO 140 1=1, MM 



flL^nrcrrJ-r) ™ 

HR=H( I , J+l) 

HOT TH (T,'J-l )+m TV3 }T*'0Y5 " 

HRH= ( H ( I » J+l } +H ( I » J ) ) *0 • 5 

HTR=TH rr+1 TJT+H (' I DT) *~0T5 

HTEMP = H{ I — 1 * J ) , 

I FlT. E'QTn RTEMP'=HT 17 JT ' 

H BH=(HT E MP +H t I * J ) )*0 # 5 
+Fl J - J hH ) ' n rvi 2 0 ,T3'0 
POINTS AT THE LEFT SIDE OF THE STEP 

TIT T F TT -TRT T 2'0TI 1 5 JT 16 

points on the bottom edge 

TrrTFTfrsTrrTE^r"r r Grrnrrrfi- 

735 HL=HL+STEPH 

H R"=H RTSTEP'R' ~ ‘ * 

hlh=hlh+steph 

■ " HRT^HRH+S'TFP’Fr ' “ * ~ 

I F C I.EQ.IH)HTH=HTH+STEP 

” TFT I”. EqTTHRTRB H=FB‘H+STE P~ ~ ~ 

GO TO 130 

ITS ‘ IT'H- I'HHTTTT » 73 5T13'0 

POINT'S IN THE POCKET 

1T7 ‘HL=HL> STEP — 

HR=HR+STEP 

‘ HCR - H CFFF5TEP" ” ~ ” 

HRH=HRH+STEP 

HTH=’HTHTSTEF r ~ ~ * " 

HBH=HBH+STEP 

GU TO‘ T'3‘0 ' 

120 IF { H — IH > 13 0 » 1 25 » 126 

"BOTTOW 'OETC) P' “C OTCfTETT 

125 HL=HL+STEPH 

HLR-HCFf+STEPH 

IF(I.EQ.IH) HTH=HTH+STEPH 
" I FT rvEQVTHRT" HBH-HBH-FSTEPH' 

GO TO 130 

I2'6 TTTI -TFfHT T 2 T *'rZ57T3 0 

points along the vertical edge 

T27 RL"=HL7FSTEP 

hlh=hlh+step 

' H TR = FTTTTT STEFT "" 

hbh=hbh+steph 

"”GO~rcr~r3'o * — 

130 IF(NDIAG) 132*132,131 

T31 ' W R'T T E( TTW 725T 

WRITE! NW,20 )HL,HR»HLH,HRH,HTH,HBH 
"'GE N E RA T ETAITITJ ) * ETC " 

132 RTEM=R( 1-1) 

*” - D‘RT= D RIT-T l — 

I F ( I ,EQ • 1 ) DRT =DR ( I ) 

- ’ iFCTTEUrr 1TRTEM=RTD 

TEMPl=Pi-AM*RU J/D22 

TEMPT = 1TO 7'PFi 0“*T DR ( I ) +DRT 

Al( I >J)=TEMP1*HL 
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ATCH O' r=TEMPT«RR" “ 

A3 (I »J)=DZ3*HLH**3/R( I ) 

' ■AV(T**J*)=az^*'HRHW3-7Rni ‘ 

A5U»J)«TEMP2*HBH**3/(DRT /RO) *(R(I)+RTEM ) 

A6TI *OT=T E MFZ* FTTH 3 7TD~R ( D7 RG T *TR" C I FFRTI + IT)' ' 

A8( I *J)=PlAM*P ( I )*CHDOT 

rF - rN - DT - A - GT '' 1 4'0~, r^'O'Tl 33 — — — 

133 WR I TE ( NW» 26 } 

WRTTETNWi 2~0TAirr;U)'»7\~21 f, J Yi A3 ( I ♦ J ) fA4 ( TVJ1 » A3TI * JT*‘A6 f r* J ) ♦ 

1 A8( I » J) 

T40 CO'NTTNUF ’ * ~ ' r ' 

KOUNT=l 

T41 DO~T4~3'~prur — — — — — - 

G ( I * 1 ) =0 * 

DCTT42 ' K=F,M‘ " ~ “ 

142 E(I»K»1)~0* 

‘ YG3~ 'CO'NTTNUE” " “ 

DO 370 J=1»N 

DO~'3TO”’ 1 = 1' * M 

I F { J.EQ.l.OR. J.EQoN) 60 TO 210 

I F C I “.“EQY 1) GO “TO 725“ ‘ - 

I F ( I • EQ . M) GO TO 210 

~~ GO"TO‘ 72'6'“~~ ‘ * “ 

725 I F ( NSYM .NE.l) GO TO 210 

T2(5 sqp=fqrt{ A‘B"surrr*j7-rnT" ~ 

SQM=SQRT(ABS(Q(I»J-1) ) ) 

SQQ=SQRTTABS ( QTTYO FT) " * ----- 

B(I) =A3 1 1 s J ) +A1 ( I » J)/ (2.0*SQM) 

'Cri)=A‘4TIVJT~A‘2'nVJ)TT2.0*SQPT ~ ----- 

ASUM=A3(I * J ) +A4 ( I » J ) +A5 ( I » J ) +A6 ( I ,J) 

F ti‘)=a'3ti »'JT*'GT(TTJ-r)'T/a nu')*sQM — — 

QTT=Q 1 1 — 1 » J ) 

I F ( I ~.FQ~. 1 TATTDT "N'SYM - « EQ. 1) QTT=TJn+r*U) 

F ( I ) =F ( I ) +A5 ( I > J )*QTT -ASUM*Q( I * J )+A6( I > J)*Q( 1 + 1* J }+A4( I , J ) *0 

1 (7 >' J+T)'-A21T'»J7^S‘QP~A8TI 7 J 1*500 ~ " *“ 

F( I.)=-F( I ) 

* "* d0“'i5O"‘k=~i;r — — ~ — 

A ( I * K ) =0 , 0 

' “T'FTK.EQTI ) AT) )“= -A“STrR-A”8'*‘0 • 5 /'SQG ’ *“ ' 

A6TT=A6 ( I » J ) 

+FC -F . EQY 1" TATJDT * NSYTT EQ'. 1 )' A6'TT=~A6 ( I »’J )+A5,CI » J f 
+F ( K.* EQ *++ 1 ) A ( I ,K)=A6TT 

I~T KTFQTr- 1TAT r7K7 = A57T, ’J') ' • ----- 

150 CONTINUE 

7 F T J .TTETORHl GO' "TO' '305 ' ' 

IF ( I.Eq.IH.OR.) I.EQ.IHH) GO TO 151 

' rFTTVGT»l HTA'NlTVmrT i'l RH7 ' GO ‘ TO T5 2 ‘ ~ 

GO TO 305 

~m ■'D'MST^5TEP^0T5' — — ~ ~ — 

GO TO 155 

" T52 D'M5P=STEP — 

155 SQQ=SQR T ( ABS ( Q ( I » J ) ) ) 

TEMP-AlTI r3T*0MSP7THri , J-T) +DMSP )" - ' 

F ( I ) = F ( I )— TEMP*SQQ 



a ~( T ~ d =xn“»l T +TEM P ^ r 0 T57'S"Q"Q “ 

305 IF(NDIAG) 310,310,306 

3'0 6 VRTTE ( N W » ZT) IV J ' 

WRITE (NW, 20 )F( I ) ,B( I ) »C ( I ) ,A( I , 1-1) »A( I , I ) *A( I ,1+1) 
' ” ~ GO" T0'”'3T0 " “ “ • 

21 0 B ( I )=0_. 

F ( I ) =0 * 

“ DO "211 KiltW ' “ • 

A ( I , K )=0 • 

i f(t. eqvk ta 1 1 VK) ■= r; ~ — 

211 CONTINUE 

G0'~rc‘~3'0‘5 = "‘~" 

310 CONTINUE 

so - 320 ' r=r,"F "* — ■ — ■ - 

DO 320 K=1 »M 

'320 Q5MA ( r,X‘I = A“( T )"+BT FJTFETf T , KTJT “ ‘ 

CALL MATINV(Q6MA,M,BB,0,DET, ID) 

~ GO "Ttf“T3~4U ,~3'30T, TD' 

3 AO DO 360 1=1, M 

" ‘ G’(T » J + 1 T= 0 » “ - - - - - 

DO 360 K=1,M 

~G ('I » UT 1T=GTT» J+ lTTQSMATlYR ) * IT ( K. ) -8 nCT*GTK» J FT *' 

E ( I , K , J+l ) =-QSMA ( I »K)*C( K) 

3“6'0 CONTTNOE" 

370 CONTINUE 

DMAs 07“ 

DO 380 1=1, M 

'D'MA=WA”X I TTONTAYAB'S fSTT7R+ X FIT 
380 DQ ( I ,N ) =G ( I ,N+1 ) 

” D'0“AW"JJ'=2',N ' ' ~~ ’ 

J=N+2- J J 

' ‘ D0'“4’00”~ T= 1~» M' " — 

DUM=0. 

DO 3 -Q--(C=1TM' - - - - - - - 

390 DUM=DUM+EC 1 ,K, J )*DGH K, J ) 

'DUM-DOMTG fl >0T "" 

DMA=AMAX1 (DMA,A8S(DUM) ) 

“A 00 *DQ'( I ,U-ll =DUM 

DO 401 I = 1 , M 

DO 40T J = IYN ■ ' 

401 Q< I ,J)=Q{ I »J)+DQl I ,J) 

+rrtTD47CG')"4'(T5TT0'5T^'0'2'-~~ ~ 

402 DO 403 J= 1 » N 
WRT T ETN ft ', 2 OT CD'Q ( TTJ JTTS 1 » MT J 
DO 404 J=1,N 

W4 WRTTETN W'i 2 0 T ( Q f T i □ ) 7TTT , M ) " 

405 CONTINUE 

"CTO TO 5ET0 ' ’ 

330 WRITE(NW,28 ) J 

■*” GO" 'T'OTOW 

560 KOUNT=KOUNT+l 

“ “ I F l FOUNT; GETEKOUNT r"GO”TO‘ ‘56 1 
IF ( DMA. GT. ERROR) GO TO 141 



561 WRITE! NW, 30 TERROR, (COUNT 

IF(NPRE.EQ.I) WR I TE ( NW » 29 ) 

DO 576 J = 1 » N 
DO 576 I = 1 » M 
P( I » J ) =SQRT ( ASS ! Q ( I.J) ) ) 

576 CONTINUE 

DO 585 1=1, M 

IF(NPRE.EQ.1 ) WRITE! NW*31 ) ( P t I » J ) ♦ J = 1 *N ] 
DO 580 J=1,N 

580 QQQ ( J ) =P ( I ,J)~1. 

NN = N-1 ' 

PP C I 1=0* 

'■ DO' 581 "J=l JUN" 

581 PP C I )=PP( I } +DTH ( J ) * ( QQQ ( J ) +QQQ ( J + l ) >*0.5 
PPC I )=R( I )*C'C*PP~( I ) 

IF (ND1 )585»585»800 
8*00 WRITET(MW»31)PP{ I) 

585 CONTINUE 
FPAD=NPAD 

mm=m-i 

WLOAD=0. 

DO 589 1=1, MM 

589 WLOAD=WlOAD+DR ( I)*(PP(1 J+PP( 1+1) )*0.5 
WUOAD=WuOAD*RO*PA*FPAD 
IFtNSYM.EQ.l )Wl_OAD = 2.0~*WLOAD 
WBAR=WLOAD/ ( 3„1416*!R0**2-RI**2)*PA) 

IF! NSYM.EQ. I )WBAR=2.*WBAR 
DO 600 I = 1 , M 

IF! I.EQ.IH.OR, T.GT. IH) GO TO 593 ‘ 

591 FRC( I ) =0. 

NN=N— 1 

DO 592 J= 1 , NN _ 

592 FRC! I)=FRC( I )+DTHIJ)7H( I ,J) 

FRC! I ) =FRC ( I )*CC*RiI )**3 

GO TO 600 

593 IF (I.GT.IHH) GO TO 591 
FRC ( I ) =”01 

NN=N-1 

DO 594 J= 1 » NN 

IF! J.EQ, JHH.OR.J.GT. JHH) TEMP=1.0/Ht I , J ) 
IF! J.LE.JHH)TEMP=1.0/!H( I , Jl+STEP) ' 

594 FRC! I )=FRC( I )-DTH( J)*TEMP 
FRC! I )=FRC( I )#'CC*R( I')**3 

600 CONTINUE 

MM=M-1 

HP0W=0. 

DO 605 1=1, MM 
605 HPOW=HPOW+FRC( I )*DR( I ) 

HPOW*HPOW*VrS*S.*28 J2TFRPS*-0**3 

HPOW=HPOW*FPAD*RPS*60. /63000.0 

HPOW=TTPOW7HMIN 

IF ( NSYM.EQ. 1 ) HP0W=2.0*HP0W 

WR I TE { NW, 34 rWL'OAO 

WRITE! NW, 35) W8AR 
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I WR~ITE CN'W >36 )HPOW 

1000 CONTINUE 

I F ( LAST )90»9Q»10 OT 
llOOl STOP 
" END 


SUB-OUT+NE MAT INVUATN1 »BTM1 »DETE- >ID) 

MATRIX INVERSION WITH ACCOMPANYING SOLUTION OF LINEAR EQUATIONS 
NOVEMBER 1962 S GOOD DAVID TAYLOR MODEL BASIN AM MAT 1 

DIMENSION A(17.17)»B(17»1> »IND£X{17»3) 

GENERAL FORM OF DIMENSION STATEMENT 

EQUIVALENCE (IROW*JROW)» (ICOLU ,JCOLU )» (AMAXo T, SWAP) 

initialization 

M=M1 
N=N 1 

10 DETER *1.0 
15 DO 20 J = 1 > N 
20 I NDEX ( J » 3 } = 0 
30 DO 550 I = 1 » N 

SEARCH FOR PIVOT ELEMENT 

40 amax=o;o 
45 DO 105 J = 1 » N 

I F ( INDEX { J » 3 ) -1 ) 60, 105» 60 
60 DO 100 KL=1»N 

+ F ( INDEX { iC» 3 ) -1 ) 80, 100* 716 
80 IF ( AMAX -ABS (A(J,K)}) 85 » 100, 100 

5 +-0W3J 
0 +COLU 

AMAX=ABS UA ( J » K ) I 

100 continue 

105 CONTINUE 

INDEXI+COLU ,3) = INDEXI+COLU ,3) +1 
260 INDEX! I'»T) =TROW 
270 I NDEX { I , 2 ) = ICOLU 


MAT 10002 
MAT 10003 
MAT 10004 

MAT 10006 
MAT 10007 
MAT 10009 

MAT 1 001 2 
MAT 10013 
MAT 10014 
MAT 10015 
MAT 10016 

MAT 10018 
MAT 10019 
MAT 10020 
MAT 10021 
MAT 10022 
MAT 10023 
MAT 10024 
MA 1 1002 5 
MAT 10026 
MAT 10027 
MAT 10028 

MAT 10030 


MAT 10033 
MAT 10034 

MAT 10036 


INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL 

130 IF (IROW-ICQLU ) 140, 310, 140 
T4~0 DETER"= -DETER ' ‘ 

150 DO 200 l=1 »N 

160 SWA'P ='AT TROW", L ] " ‘ ” 

170 A ( I ROW »i_)~A( ICOLU »L) 

200 ATICOru , L) =SWAP 

+F(M) 310, 310, 210 


MAT 10038 
MAT 10039 
MAT 10040 


MAT 10043 
■*iA T 10044 


T 100<+7 





- 98 . 


~2T0 DO ~ 2'50 'T=T» M~ 

220 SWAP=B(+-OW,L) 

230 B t I ROW * i_ ) = B ( I COLU »L)' 

250 B ( I COLU *L)=SWAP 

DIV+DE PIVOT ROW BY PIVOT ELEMENT 

310 PIVOT =A(ICOlU *IC0LU ) 

— - D ETE r -=d*e r e r>p r v ot * 

330 A ( +COLU » +COLU ) = 1 • 0 
"3 VO' DO" 33'0" "L= 1 » N ' ” ' ~ * 

350 AIICOLU > L ) =A ( I COLU *L)/PIVOT 

353 +FKM) '3F0*T~'3'8Q 3"67) 

360 DO 370 l=1»M 

"370 B frCOLU » L~) =B ( I~COLU »L)7 PIVOT 

-EDUCE "NON— P I VOT -OWS" ’ ' 

~3*80 D'O "550' L1'=1',N 

390 IFILI-ICOlU ) 400, 550, 400 
'400 T=A“( LI ,'ICOLU ) 

420 A ( LI > I COLU ) =0 • 0 
' I Ft T T430,55'0 ,^3D 
430 DO 450 i_= 1 » N 

'~V5"0 ATL1»D =A"(L1 ,l )-A( ICOtU’ »L )"*T 
455 I F ( M ) 550, 550, 460 
" V60 DO 300 U=1 ,*M 
500 B(L1»L)=B(L1,L)-B( ICOLU »L)*T 
■ 55"0 'CONTINUE * 

rNTER-CFTANGE- 'COLUMN'S " 


' 600 D'0'710 T=1,'N 
610 L=N+1-I 

~6‘2'0 IF ( rND'EX(Lil')-I'N'DcXlX,2>') "630, 710T630 
630 JR0W= INDEX(L,1) 

" 6"4'0 J CO LU” =TfTD EX T L , 2 T “ 

650 DO 705 K= 1 , N 
' '66>0 SWA'P=A"0C, JR'OW) 

670 A ( K , JRO W ) = A l K * JCOLU ) 

70*0 A (X, JC’OlU )'=SWAP 
705 CONTINUE 

~7T0 CONTINUE* * 

DO 730 K. = 1 »N 

I*F ( I N DE’X't (7, *3 )' -T) '7157720,715 
720 CONTINUE 
“7'3"0 'CO'NTI'NUE' ----- 
ID=1 

T40 -ETLRN ” ----- ~ 

715 ID =2 

~G0' TCT’74'0' 

END 


MAT 10048 
MAT 10049 


MAT 10 05 2 
MAT 10000 
M'AT10053 


M'AT 10057 

MAT 10059 
MAT 10060 

MAT 10062 
MAT 10063 
MAT 10064 
M'AT 100*65 


MAT 10069 

MAT 1 007 1 
MAT 10072 

MAT 10074 
MAT 10075 
MAT1007'6 " 
MAT 10077 
MAT 1 0078 
MAT 10079 
MAT 10080 
MAT 10081 

MAT 10083 
MAT 10084 


MAT 10087 
'MAT 10088 
MAT 10089 
MAT 10090 
MAT 10093 
MAT 10094 

MAT T'0096~ 
MAT 10091 
MAT 100^2 * 



Card 1 Format (72 H) 

Identification Card 
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Card 2 Format (1615) 

Last - Integer to determine whether additional input data 
are to be read 

Last =1 , no more input data 

Last =0 , more input data to be read from state- 

ment 90. 

NPAD - Number of pads (see Fig.A2) 

M •• Nunber of grids in the radial direction 

N - Number of grids in the circumferential direction 

IH Grid number for the bottom edge of the step (see Fig. A2) 

For NSYM=1 } set IH-1 

IHH - Grid number for the top edge of the step (see Fig. A2) 

JH - Set JH*1 

JHH - Grid number for the left edge of the step 

LKOUNT - Maximum number of iterations allowed for the pressure 

calculation (recommended value: 10-20) 

NDIAG - Control for diagnostics 

NDIAG =1 , diagnostics print out 

NDIAG =0 , no diagnostics 

IRRG - Control for irregular grids 

IRRG =1 , read in irregular grid spacings 

IRRG =0 , uniform grid spacing 

NFRE - Control for printing out pressure profile 

NPRE = 1 , print out pressure 

NFRE =0 , no pressure print out 

NYSM - This integer is used to control whether pressure 

calculation is made for a full pad or half a pad as 
in the case where the pressure is symetrical about 
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the center line (see Fig,A2). 

NYSM =0 , For calculation covering full pad 

NYSM a i , For symmetrical pressure profiles where 

where calculation is only made for half a pad. 

ND1 - Set ND1 =0 , for normal runs. 

Card 3 Format (8F10.7) 


DO 


Outside diameter of thrust brg, , 

d in Fig.AI (in.) 
o 

DX 

- 

d. in Fig.AI (in.) 

X 

THG 

- 

Angle extending the pocket region, 6 in Fig.A2, (degrees) 

S 

THL 

- 

Angle extending the land region, 0^ *- n Fig.A2, (degrees) 

STEPD 


Depth of the step (in.) 

WO 

- 

Outerwidth of the shreud, W q in Fig.A2, (in.) 

WI 

- 

Innerwidth of the shroud, in Fig.A2 (in.) 

ERROR 

- 

Convergence factor for pressure iteration, 
(recommended value: .001 - .0002) 


Card 4 Format (8F 10,7) 

This card is required only when 1RRG = 1 

DR(I), I a 1, MM • Dimensionless irregular gird spacings in the radial 

direction. (MM - M-l) 



Fig.AI 
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Card 5 


Format (8F 10.7) 


This card is required only when IKRG = 1 


DTH(J), J = 1, NN - Dimensionless grid spacings in the circumferential 

direction, (NN = N-l) 


(IHH,JH) 



UH.J H )V^ 7yr~^0 H) f x 

(•JluiUTt,*') Fig. A- 2 




DTH(J) = 


(A0), 


0 L + 0 C 


9 J + i ~ 8 i 
®L + ®G 


Card 6 Format (15, 5X, 7E 10.3) 

NVTSM - Total number of viscosities to be investigated in the 

production run. 

VISA(I), 1=1, NVISM 

2 

The arrary of viscosities, (lb-sec/in ) 

Card 7 Format (15, 5X, 7E 10.3) 

NKPSM - Total number of angular speeds to be investigated in 

the production run. 

KPSA(I), I - 1, NRPSM 

The array of angular speeds, (Rev. per sec.) 

Card 8 Format (15, 5X, 7E 10,3) 

NPAM - Total nunber of ambient pressure to be investigated 

in the production run. 

PAA(I), 1=1, NPAM 

The arrary of ambient pressure, (PSI) 
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Card 9 Fermat (15, 5X, 7E 10.3) 

This card reads in the film thickness to be investigated in the 
production run. In the case of a parallel film, the film thickness 
in the land region will be read in the array designated as HMINA(I). 

In the case of non-parallel film with an axisymmetric coning of dishirg, 
BMIKA (I) represents the film thickness at the inside diameter in the 
land region. The film thickness at the outside diameter is HOUT(I) 
xdiich read in Card 10. 

l 

l 

The variables to be read in this card are: 

HHMM - Total number of film thickness to be investigated in 

the production run. 

BMLNA(I), 1=1, NAMM 

HMINA is the film thickness at the inside diameter in 
the land region (inch) 

Card 10 Format (8F 10.7) 

H0DT(I), I = 1, NHMM 

HOUT is the film thickness at the outside diameter in 
the land region, (inch) 

Card 11 Format (15, 5X, 7E 10.4/(8E10.4) 

NHDTM - Total number of velocity or time variation of gas 

film thickness to be investigated in the production run. 

HDOT(I) - Array of velocity or time variation of gas film thickness 

in in. /sec. 
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P R 0G RAT 5 !" R'S'G'AO (* I'N POTTO U'T PUT)' ” ' 

DIMENSION WN ( 30 ) » QN ( 7 ) , AG(9)»A0G<9) 

D+MENS+O'N AHAtlOOD* FXaU'OO)* F2A‘(10'CV» F3*A1 lOOTi ~F4 A ( 100) » FII9) 
DIMENSION F5A ( 100 ) » F6 A ( 100 ) » F7A { 100 ) * F8A(100)» F9A(100) 

THE TATTuES ‘ARE PREPARED"" I N THE FOLLOWING ORDER THAI (TJ'AII ) >J=1T9) 

= ( Dl*ll» D2*I1A> Dl* I 4 > D2*I4A» -D2*I1A0» -D2*I+AO» D2* 1 5/H0**2 • 5 
X"DI'*r5A7HOW2.5X ^D'2TfT5A0VH0*Tf-2;'r'I'T“ 

DATA F1AI 1 }/0./»F2A( 1 ) /7 . 854/ * F3A { 1 ) / 3 . 1416/ * F4A < 1) /0./ 

DATA CFTA ( I P* I = 2T60 YT 7.8 5552 8 1 £- 02 XI* .“5 720 3 4 XE^Ol X 2. 36 03 750 E-0 1 » ' 
2 3*15151 2 3E-0 1 » 3. 9463912E-01 * 4.7459680E-01 ♦ 5 . 5512133E-0 1 » 

2 *6 .36 311 4*6 E-0 IV “7.1 8258015= 01 X 8T0TCT941T5-01 »”8X8'4"8*956*0E=01 , ' 

2 9.6978133E-01 » 1 . 0558636E+00 > 1 . 1432584E+00 * 1. 232086 1E+00 » 

2 ' 1 .'3 224*7155+00+ “1 ;4X45“448 E+ CTOT ' 1 .*50344 1 5E+0'0+'**T.5'T)A3'034*5+ 0 0 X~~ 

2 1.7022791E+00* 1.8025243E+00 » 1 .905203 1E+00 , 2.0104879E+00 » 

2 ■*2+1185609E+00+~ 2;X2*9"6*l*45E+*00» 2 . 3 435 52*9 5+0 0 XX . 4 61 49 1 9 E+ 00 ♦ 

2 611E+00I 2C707 04 E+OOT 2C837183 5E+00 » 2. 708745E+00* 

2 3 i 10 2 7 4'8 E+Q 0 X 3+2 5X7 020 E+0 0 » 3.40X49 64 E+0 0+ “3 .5 5'6 Q2 3 76 + 0 0 * ' 

2 3 . 7166766E+00 1> 3C 8387 1E+00* 4C0580843E+00 » 4.23 7876E+00* 

2' 4+42 '52*0*65+0*0 X“4C6ZT85025+*O0"f 4 . "3543 2 15+ 00X5+ 05-29163 5+00 . ” 

2 5. 2810529E+C0* 5«5206483E + 00» 5 . 7725837E+00 » 6 .0378224E+00 * 

2 6+ 3 1*7420 3 E*+X 0X5 « 6X25372 E+00+ 6.9 24447 65+0 0 T~7 .25455706+00+ ' 

2 7 .6044175E+00 * 7.9757465E+00 » 8 . 370449 1E+0Q » 8 . 7906429E+00 » 

2 9". 23*85'8 745 +‘00 X 9„7 172 1‘83 E+00 , 1.0229187E +011TT. 077790 85+ 01/ 

DATAUF1AU+) *1=61*100}/ 

2 ’ ■ 1 .T3 67IX4E+01 V 1+ 2*00*1*0 1 85+0 1 i 1.2 68*43 9 6E+0 1 * 1 • 3422 6 7 35+ OT » 

2 1.4222038E+01* 1„ 5089573E+01 » 1 , 603341 2E+0 1 * 1 . 7062939E+01 ♦ 

2 1 . 8T 89023 E+’Ol ?* 1 . 942430* *£+0 1 » 2 . 0 7 8 3 5 8 0 E +0 IX X . 2264 2 0*65+ 01*" 

2 2.3946707E+01» 2 • 5795465E+U 1 » 2 . 7859629E+0 1 » 3. 0174285E+01 » 

2 *3*.278'1955E+0TY‘3.5T346276+0T» 3 . 9096 2 8 7 E+ 01 » “4 . 2 746 5 3 06+01 * 

2 4.7385298E+01 , 5 • 253947 5 E+U 1 » 5 . 857209 5E+0 1 * 6. 5695462E+01 * 

2 “ ' 7X4 19*0157 E+OT P'S*. *44332 1 25+0 1 » 9 .‘6 94*08 26EX0TXT .1X434 925+ 02 X"' * '*" 

2 1.31 5000E+Q2, 1.57010 5E+02, 1 . 899410 1E+02 * 2.3442157E+02 , 

2 * 2 • 9656997E+02 X 3+ 87X653* E+02 * 5. 2*6*667 4*45+02‘X7**. 5*7 8*89 7 6E+02 X 
2 1.1832960E+03T 2. 1018302E+03 » 4 . 7245976E+03 , 1 . 887838 1E+04/ 

DA T A IF 2 A* (1 f»3 =2+5*0 T7~ 7 .'8 58*6 2 1 5 E +0 0 X *7 . 8 725 5 5TE+ 0 0 X“7 .*8558 2 3*95+0 0 * 

2 7.9284975E+00* 7 .9706737E+00 , 8 .0224795E+00 , 8.08,407125+00* 

2* ”8 ".'1 55556*05+ 0737" '8+23*7 3'9 22 E+ 0 OV * 8 3X95 9T45+0 CT~8 . 4 3X 51935+ 0 0 i ' 

2 8 • 546 497 6 E+00 » 8 . 67 1 8 859E+00 * 8 . 8090839E+00 » 8 . 9585338E+00 * 

2 '9. 12*0723*05+00 X 9+ 2*9*61*87 15+00 » 9 . 48551 555+ 00X 9 .53934 515+ 0 0 X" 

2 9.9083843E+00* 1.01433 8E+01, 1 .0395223E+01 * 1 . 066477 2E+0 1 » 

2 1. 09 530 375+01 » 'T;iX6Tr0'4E+0*l> 1 . 1*5*90 15 25+OXV' 1X194X4675+0 IX ~ 

2 1. 2316453E+01 » 1 . 27 16640E+0 1 » 1 . 3 143699E+0 1 > 1 .3599453E+01 » 

’ ' " 1C40 8 5“ • TE+’OTTl (7450 5 2075+0 I ; 1C5X59775E+OTX I .*575' Z2T” 5+ 0 T * ***' 

2 1 . 6385407E+0 1 » 1 .7062493E+01 » 1 . 7786946E+0 1 » 1 . 8562586E+0 1 » 

2 ~ 1*. 935*35245+01 ♦“"ZYUZB+TTl 2'E+OT , * 2 *.12*40*9 95E+0 1 * ~2 . 2258 1705+0 1 i * 

2 2.3372563E+01 » 2.4561200E+01 , 2 . 58419 10E+0 1 » 2.7223422E+01 ♦- 

2 — 2 ."87T5497 E+ OT » 3 ."032 9 '0*686+0 1. 3*. 237*64*0 5E+*0 1 X 3.35Tr315F+0TX ' * 

2 3.602 365E+01 » 3 . 826 8 1 50E+U 1 ♦ 4 . 07076 1 3E+01 » 4. 3370409E+01 » 

2 4 .5282353 E'+ 01 »~4 "."947 2*9*4 l'E+0 1 ♦ 5 . 2 9759835+ 017 “5*. 683 03 3*9 E+ 01 / * * “ 

DATA{F2A( I ) »I=61T100)/ 

2 6* 1080'85'6'E+Cl » 6 . 5T79435E+01 » 7 . 0 9 8 6 3 7 9 5+ 0 1+* 7.6 772 0 24*E + 0 1 

2 8.3218737 E+ 01* 9.04Z3387E+01 * . 850041 6E+01 * 1 . 0758568E+02 * 

2 *l‘i X784127E+02 » ” 1 .29461 60E+U2 » ”1 .4X6'8123E+0 2» X . 5778484E+02 » 
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2 1.75 1202 1E-RJ2 , ~T.95T 1577 E+02 , 2 .16299 9 AE+UZ -"2.4 5776 90E+ 0 2 » 

2 2.77 0607 1 E+02 * 3.14532 6E+02* 3.5911751E+02 » 4. 1258599E+02 ♦ 

2 49 772'6736E + Q2 * 5 .5‘6'2'6 38 2E+02 * 6 . 53T6773E+02 » 7. 7553624E+Q2 » 

2 9.29617 5 5 E+02 > 1.12748 5E+03, 1 . 38589 38E+03 > 1 . 72988 27E+03 . 

2 2 • 19506 10E+03 > '2.85I9737E+03 » 3 . 7937283E+03 » ‘5 • 2009564E+Q3 * 

2 7C400 326E+03, 1C1 04100 1 E+04T 1C7522447E+04 » 3.0261024E+04, 

2 5.90686" 8'4"E+“0 4‘» T .’3 993T3'5'E+ 05 , 4 . 7 198425E9-~05 , "3. 7735728E+06 / ' 

++ OF A0++0+ 2060D / 3 . 14195 10E+00 > 3. 1429824E+00 » 3. 1447027E+00 » 
2 3 * 1471 14ZE+00 * 3 . 1 50'2 198E+00 » 3 . 154D23~5E+00 ♦ 37 1585298E+00 . 

2 3 • 16 37446 E+00 > 3 . 1696743E+00 » 3 . 176326 5E+00 » 3. 1837096E+00* 

2 3.197 873 TE+Q 0 » ' 3 .‘ 2*00707 5 E+UG » 3 . 2T03442E+00 »' 3 . 22075 58E+00 » 

2 3. 2319562E+00* 3 .24396 02E+UQ » 3 . 2567839E+00 » 3 • 2 704448E+Q0 » 

2 3 • 2849616 E-FQO j 3 . 30W5"46E+00 » 3 . 3T6*6'4 5 3 E+ 00 77 '.'33385 6 9 E+ Q 0 » 

2 3.352U144E+00 » 3 . 37 1 1442E+00 » 5 . 3912749E+00 , 3 . 4124369E+00 ♦ 

2 3 +43466 2 5 E+0 0 * 3 .4579865E+00 » 3 . 4824-45 8ET+00 » "3 . 5080801E+00 ♦ 

2 3.5349315E+00 » 3 . 5630451 E+00 » 3 . 5 92469 1E+00 » 3 . 6232549E+00 » 

2 3 * 6654576E+00 » 3.6891360E+00 » 3 .7243531E+0Q > '3. 7611763E+00 » ' 

2 3»7996778E+OOj 3 .83993 52E+0Q , 3.8820315E+00 * 3.9260561E+00 » 

2 3 « 9 72 1049 E+00 » 4'.02U'2 _ 8l2E+00 . 4 . 07W963E+00 » ”4 . T274-7Q2E+00 ♦ 

2 4. 1787323E+00 , 4 . 23662 28 E+00 , 4. 2972932E+00 » 4. 3609078E+00 » 

2 "4.4276449E+0 0 t 4i497“69 84E+00 » 4 . 5 712794E+0Q > '4 . 648*6 182E+00 ♦ 

2 4 . 7299663 1+0 0 » 4.81559 1E+00* 4 . 9058 1 8 8E+00 » 5.0009572E+00/ 

DAT A ( F 3A' ( I ) » I =6 1 » TOOT 7 * 

2 5. 10138U1E+00* 5 .207491 1E+00 » 3 . 3 197370E+00 » 5 . 43861 34E+00 * 

2 ' 5.5646716E+0D7 57 6*98 5Z67 E+ 0 0 i 5 . 84Q3664E+007 "5 . 9924*634*E+00 » 

2 6.1541 7 5 E+C 0 s 6.32702156+00 » b . 5 12080 6E+Q0 » 6 . 7 106352E+00 » 

2 6 • 9241382E+00 » 7. 174-Z599E+00 * 7 . 4029288E+00 s» "7.6723848E+00» 

2 7.9652429E+00, 8 . 2845 764E+00 » 8.6340204E+00 » 9 • 0 179069E+00 » 

2 .44143 0 E+00 > 79T09212E+0Q > 1 . 0434067E+0 1 > 1 • 102*0413E+0 1 » 

2 1.168 1892E+01 » 1 , 2433635E+0 1 ♦ 1 . 3295 1 14E+0 1 * 1 . 42918 18E+0 1 » 

2 T.5457768E+01,- I7683"9404E+01, 1 . 85U1808E+0 1 ♦ 2. 0539089^+0 1 ♦ 

2 2. 3092560E+01 , 2 .63845 28E+0 1 , 3 . 07859 14E+0 1 » 3 . 6965 1 37E+0 1 » 

2 4.6260583E3-01 , '6„1T9‘8892E+01 j 9 . 29T1744E+0 1 » 1. 8680213E+02/ 

+AT1UF4AU+DT+32T60D/ 6. 736515E-02, 1.3755 14E-01» 2 . 0655024E-0 1 » 
2 2.75795718-01, 3 ,.45'3902TE-01 , 4. 1541249E-0 1 » 4. 8595455E-0 1 . 

2 5.57107 16 E— 0 1 > 6 » 2896280E-0 1 » 7 . 0 16162 5E-0 1 ♦ 7 • 75 1648 1E-0 1 » 

2 8. 497086 3E-01, "9. 2535 1 02E-01, 1 . 002198 9E+00 8 1 . 080363TE+00 » 

2 1 . 1599589E+00 » 1 .241 1 062 1+00 » 1.3239304E+0 0, 1 . 408 5626E+00 * 

2 7 .'495140 1E4-00 s''T.583'8Q 73E+00 > 1 . 6747 16 1E+00 » 1 . 768 0265E+00 ♦ 

2 1.8639077E+00. 1 .9625385E+00 » 2 . 064108 5E+00 » 2 . 1688 189 E+00 ♦ 

2 2 • 2768835E+007 * 273885300E+O0 » 2.504D014E+003 2 • 62755707+00 > * 

2 2.747 47 42E+00» 2.8760502E+U0 , 3 . 009603 8E+00 * 3. 1484776E+00 » 

2" 3;2'97040Q"r+70r71447'6~882E+00* 3. 6008 5 07 E+ 007 3.77499G7E+00 » ” 

2 3.9366095E+0Q, 4 . 1 16 25 04E+00 j 4, 304503 7E+00 » 4 . 5020 1 10E+00 » 

2 4.70'94711E+00» ~4.9‘2T6467E+00, 5 . 1 5737 1 3E+00 » 5 . 3795 573'E+OO » 

2 5.6552054E+O0i 5 . 92541 54L+00 , o . 2 1 lb97 8E+00 , 6 . 5 144879E+00 ♦ 

2 "6.8 3~616 20*E+”0 0 » 7 . 1T80548E+00 , 7 . 547981 2E+00 » 7 * 9299605E+00 » 

2 8 . 3442447E+00 > 8.7873513E+O0 » 9 . 262 102+E+00 > 9 . 7716694E+00 / 

D A rA7T4ATTT7"I = 67 71 071 ) / 

2 1C031 627E+01, 1.0 10015E+U1T 1C1547415E+01 * 1 , 2237038E+01 . 

2 1 *‘29B"48 3ZE+0 1 » 1.3797611E+01 , 1 ,46832 10E+01 » 1 . 5650673E+01 »” 

2 1 • 67 10487E+01 , 1.7874863E+01 , 1 . 91 5809 1E+01 » 2 . 0576979E+0 1 » 

2 2.*2T5139*8'E+01 » 2.390498T6E+01 , 2 .5886029E+01 ♦ " 2 , 80'68606E+0 1 » 

2 3C0554O7OE+O1, 3C3372 1E+01T 3 C658773 lE+Ol * 4.0275 05E+01, 
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2 4.4535078E+0X 

2 7.0377289E+01 

2 1.2647590E+02 

2 2.8763837E+02 

2 1.1634881E+03 

DATA ( F 5 A C I ) ,1 = 1 
2 8 . 26968463-0 X 

2 1.9536883E+00 

2 3.13 8161 E+00 

2 4.4239953E+00 

2 5. 85 17341 E+00 

2 7.4791 17 1E+00 

2 "9.37 797 57 E+0 0 
2 1C1643017E+01 

2 1.44033367+01 

2 1C7 38 77E+01 

2 2 • 2207604E+Q1 

2 2 • 788 9564E+0 1 

2 3 . 546 1763 E+01 

2 4.5830332E+01 

DATA ( F5'A ( I J ,1 = 60 
2 6.047467 1E+0 1 

2 8 . 1923058E+01 

2 1. 1473976E+02 

2 1.6773775E+02 

2 2 . 5 9438 3 8E+02 

2 4 • 3‘3'3 0 9 2 6 E +07 

2 8.0793374E+02 

2 li78'6085 2E+03 

2 5 • 3189242E+03 

2 3 » 06426 97 E"+0"4 

2 3.783U331E+06/ 

DATA1P67U n , r=r 
2 7 .86795 3 3 E+00 

2 7 • 9303414E+0Q 

2 8 . 044542 2 E+00 

2 B.2T39'440E+0'0 

2 8 • 443738 9E+00 

2 8C74I2 3E+00" 

2 9.1 166732E+Q0 

2 .583-502 E+00 

2 1.0 160 104E+0 1 

2 1 • 0 87T2 8 7E+ 0 1 

2 1 . 17 51005E+01 

Z 1.2 ' 46"4‘67E+01 
2 1.4224658E+01 

2 1 .598702"4E+aT 

DATA ( F6A( I ) ,1 = 60 
2 1.826 76 783+ 'TT 

2 2.1305867E+01 

2 -2.5-4-671 09 E+OT 

2 3 • 1386243E+01 

2 4C02324-17E+01 

2 5 .43627 33E+01 


4. 94892 38E+0 1 
8*0293107 1+0 1 
1.5091 947E+02 
3.7671568E+U2 
2.0745753E+03 
» 59)/ 0. 
"1.T050'180E'+00 
2 <t 2432295E+00 
3.45W33C3+00 
4. 765 49 45 E+00 
6. 23737457+00 
7.925131 2E+00 
9'. 995 50 9 7 E + 0 0 
1C2280673E+01T 
1 . 5T90438E+0 1 
1C 831251E+01T 
Z* 34865 60 E+01 
2.9576990E+01 
'3. 7 7 46313 E+ 01 
4.9014836E+01 
, 100 )/ 

6 • 50662 98E+0 1 
- 8.88160-87E+01 
1 .2560462E+02 
'1 » 8"5"95~3 37E+0Z 
2.92521 05E+02 
5 . 0027714E+02 
9.6635414E+02 
■ 2 V2641687E+03 
7.5502829E+03 
- 5 - . 9616 41 7 7 E + 0 A , 


i 59)/ 7.654018 
7. 8788175 E+00, 
7. 95 3 93 00 E+00 , 
8C0815478E+00, 
■'8'. 265 48 3' E+00, 
8.5 1 1432 1E+00 , 
C ‘27 3737E+OOT 
9.2241836E+00, 
C716 52 053+00, 
1.0324138E + Q1 , 
"7 • F07 38383+0 1 , 
1.2002430E+U1, 
IC3T5l*32 "E+01 , 
1.46239 16E + 01 , 
“T. 64976 7E+01 , 
,100)/ 

~r. 89^62 233+ 01 , 
2 . 222 197 1E+01 , 
- 2 . 6-7'4- 86CT0E + 0 1 , 
3.3259352E+01 , 
9-C31342 01E+01T 
5 . 923 1660E+0 1 , 


5. 529787 2E+01, ~6. 
9.2423422E+01, 1. 
1.83X012QE+02," 2. 
5.1416995 E+ 02, 7. 
4.6 82 2 3 3 3 E+ 0 3', 1. 

, 2.7497442E- 
1 . 3851 3 1 57+ 0 0",~ 1. 
<£•5370353 E+00, 2. 
3C7670041E+00 »" "4. 
5.1 167 138E+00 , 5. 
6.6364245’E+OOr 7. 
8.3893125 E+00, 8. 
1.04575887+013 1. 
1C295 147 3E+0 IT 1. 
1.6072 81 8E+01, 1. 
1C 886382E+0 IT 2. 
2.4854161E+01," 2. 
3.1392586E+01, 3. 
4. 022 14663+01* 4. 
5 . 2492622E+0 1 , 5. 


2 r68 3737+01 * 
0747819E + 02 , 
2665762E+02, 
4249041E+02, 

8 78 9 9 67 E+ 04/ 
01,5. 5045962 E-01 
667840"8E+00 » ' 
8356944E+00 , 
0914178E+00, 
4784962E+00 , 
0499568E+00’, 
8730716E+00, 
r03"6066E+01 , " 
3658065E+01, 
69042507+0 1 , 
1009798E+01 » 
6318738E+01, 
3349322E+0 1 , 
29-08234E+01, 
6299036E+01 / 


7 • 0 128000E+0 1 , 

9 . 6499 133E+0 1 »~ 
1.3788928E+02, 
2. 069 2 2 747+ 02 T" 
3.3152694E+02, 
5.8192238E+02,- 
1 • 169449 1E+03 , 

2 * 9307 839E+03 , 
1.1236044E+04, 

1 .409-879 8E+ 05, ‘ 


7 . 5 7 22 5 44 E+01 , 
1. 0509380E+02 , 
1.5183581E+02, 
2. 31T91Q3E+Q2 » 
3. 7785238E+02 , 
6.8251978E+02, 
1.4342207E+03, 
3.8 893 080E+03, 
1.7787769 E+ 04, 
4. 7 435 6 06 E+ 05 , 


4E+003 7.S555648E 
7. 8928191 E+00, 7 
7. 98079747+00 »”8 
. 1 220643E+00 , 8 
8. 3208 88 5E +"00," 8 
8 . 5 8 348 60 E+00 , 8 
C 164 5ZE+00, 

9 , 3376399E+00 , 9 
. 567 0 9 33+ 00"*" 1 
1C0497005E+01, 1 
1. 1787 50 13+0 IV 1 
1 • Z268 144E+0 1 , 1 
1C349'501 IE+0 1 » 1 
1 . 5048615E+01 , 1 
1 • 7U477887+0 1 *"' 1 


+ 00 ,7 .8602071 E+0 0 
• 9099834E+00 , 

.0 X097 7 93 + 00,“ 

. 1661682E+00, 
.3807 6726+00,' 

• 66005Q7E+0Q , 

• 0l"48 5 5 5'E+00 , 
.4573150E+00, 

• 000'4438E+0 1 , 
•0679206 E+01, 

•“1 57 2 9 703+01 , 
•2549132E+01 » 
.37489323+01, 
•5500868 E+01* 

. 7636 5667+01/ 


1 . 9673 - 83777017 " 

2.3214764E+01, 

2 .8X512097+0 1 »' 
3.5 336290E+0 1 , 
4C640 644E+Q IT" 
6 . 48635 a 6 E+0 1 , 


2.0'4587403+01 , 
2. 429 3085 E+01, 
2.9690866E+01, 
3.7648184E+01, 
5. 0723481 E+0 1 , 
7 a 142942 1E+Q 1 » 
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2 7.9 rs 119SEV01 »~8“. 83ZT633t:+01 » '9.9333 2 2'8'£+0 1' *~t . 12721 9 32+ 0 2 » 

2 1.2 23557E+02, 1C4 43216+02, !C76416b3E+02 » 2 . 1 1 0456QE+02 , 

2 2.5760616E+02, 312235870E+02 » 4. 1630690ET+02 » "5 .6O28454E+02 , * 

2 7.9777874E+02 , 1 . 2326001E+03 , 2 . 166835 3E+03 » 4. 8210180E+03 » 

2 1. 015 072E+04/ ‘ .... - 

DATAUF7AU+) »+=l* 3 )/ 6C2 32000E+00, 6 . 2845747E+0G *6 . 28870 12E+00 » 
2 6. 2'955873'E+O'O » 6 . 3 05 2457E+0 0 » 6.31769 45E~+0 CJV "6 . 3 3'2'9'57'CTE+ 0 0' » " ‘ ' 

2 6 • 35 106 18 E+00 » 6.3720431E+00 , 6 . 3959406E+00 * 6 .4227999E+00 * 

2 6.4526727E+00, 6 . 48“5'6 166E+00 , 6 . 52T6958E+00 »' 6 •5W98T5'E+00 * * 

2 6.60355 14E+00 » 6C64 4 11E+0Q, 6 ,6988937E+00» 6.7518606E+00 » 

2 6 • 8 08 5'0 1 8 E+0 0 » 6 . 8'689’366‘E+OO » 6 .93'32940E+00 7.00T7134'E+00 » 

2 7.0743452E+00. 7 . 1 5 135 19E+UQ , 7 .2329084E+0Q , 7 . 3 192034E+00 » 

2 7.4 TO 44 012+0 0 » 7 . 5‘0W37‘6E+0O , 7 . 6 07T6 3 1 92+ 0 0 » "7 i7T6 0 7 7"4E+ 00*’ 

2 7 • 82 9448 2 E+00 » 7 . 949040 1 E+00 , 8.0751721E+00, 8 . 208 188 8E+00 > 

2 8. 3484623'E+OO , 8 .49&3949E+00 » 8 .6524219E+0G »' '8 . 8 r70147E + 00 » ' 

2 8.9906842E+00* 9 . 1739847E+00 » 9 . 3675 184E+00 » 9.5719399E+00 » 

2 9 . 78 796 2 2 E+00 » 1 . 0U16362 E+01 , 1 . 0257989E+0 1 1 . 05X37692+0 1 » 

2 1. 07847 18E+01 , 1.1071 4 E+01» 1 . 1376685E+0 1 » 1 . 170027 1E+0 1 » 

2 1 . 20 44 190 E+0 1 , 1 . Z410079E+U1 » 1 . 2T9975 1E+0 1 » ' 1 .32X5214E+0 1 » 

2 1 . 36 5 87 00 E+0 1 » 1 . 4132695E+01 > 1 .4639974E+01 » 1. 5183642E+01/ 

D A T A CF 7 A ( I T»I~ 60,100)/ 

2 1 . 5767 183 E+0 1 > 1 .6394521E+01 > 1 . 7070076E+0 i » 1 . 77988 56E+0 1 * 

2 1.858-654&-E+01, 1 .9439152 1 E+0 1 » 2 . 032548'62+0 1 V 2 . 1372637E+0 1 »“ 

2 2. 247086 3 E+0 1 » 2.36714 1E+01* 2 . 4987684E+0 1 s> 2 . 64348 1 1E+0 1 » 

2 2. 030 11E+01, 2, 7 726 E+Gl, 3C1759 145E+0 i ‘3 . 39467052+0 1 > ' 

2 3.63962 iOE+Ol , 3 . 9 15 1 5 65E+01 » 4. 2266341E+0 1 » 4. 5806456E+0 1 , 

2 4C 53764 E+0 IT 5C4310 252+01 , 3. 90 1 112+0 XT" 6. 621237 32+0 1 » 

2 7 . 364097 5 E+0 1 , 8.24806926+01 , 9 . 3116396E+0 1 * 1 .0607436E+02 , 

2 1.22 08965 E+02, 1 .422143TE+02 , 1 .6T9971ZE+0Z, ' 2 . 0T79468E+02 » 

2 2 . 47 3 36 60 E+0 2 » 3 . 1 08 1 241 E+02 , 4. 03 1 146 3E+02 , 5.4489157E+02 » 

2 7C7 2 615E+02T 1C2CT 46 E+03, 2C13'59358E+03T~ 4. 7745321E+03 ' 

2 1.897567 1E+34/ 

DAT A TF8A ( I } » I = 1 » 59*)/‘0.“ ' 2.7497442E-Qr,5.30:45962E-Ql » 

2 8.26 6846E-01, 1. 1050180E+00 , 1 .3851315E+00, 1 . 667840 8E+00 » 

2 1 • 953688 3E+00 » 2.2432295E+00 » 2 .'5370353E+00 2 . 83369442+00 » 

2 3.13 8161E+00, 3 . 4500 330E+00 » 3 .7670041E+00 , 4. 09 141 78E+00 » 

2 4. 4 2399 5 3 E+ DO , 4.T634945E+00 , 5,11671382+00,“ 5 .'4T8496ZE+00 

2 5 * 85 1734 1E + 00 , 6. 2373745E+00 , 6 .6364245E+00 » 7 . 0499568E+00 , 

2 7.479117 1E + D0 , 7 . 9251 3 12E+00 , S.3833125E+00, 8 . 87307X62+00 »’ 

2 9.3779257E+00, 9 . 9055 097E+00 , 1 .0457588E+0 1* 1 . 1036068E+01 » 

2 1. 16430 17 E+0 1 , 1 . 22&0673E+01 , 1 . 2951473E+0 1 T 1.363806-52+01, 

2 1.4403336E+01, 1 . 5190438E+01 , 1 .60228 18E+01 , 1 . 6904250E+0 1 » 

2 1.7 83'8 8T72+ 01 » “ 1 . 8*83 I2'51 E+ 01 » 1 .'98'8’6382E+01 V 2. TOO 9 7 922+ 01 »" 

2 2.2207604E+01, 2C3486560E+01 , 2 .485416 1E+0 1 » 2 . 63 18738E+0 1 , 

2 2.7 839 5 6 42 +01, 2 .'937699OE+01 , 3 . 1392526E+01 »“ 3. 33493222+0 1 » 

2 3.5461763E+01 , 3 . 77463 1 3 E+0 1 , 4.0221486E+01 , 4. 2908234E+0 1 ♦ 

2 4 • 58303322+01 »'4*9'0 14836 E+Q1» 5 . 2492622E+0 1 » 5 .“6'299036“E+01 / 

DATA ( F8 A ( + ) ,1=60,100)/ 

2 'H‘.0474'67TE*+01 '6.5'0 - 6'62'9'8E+01 , 7 . 0T2 800 OE+013' 7«57225*44E+01 > 

2 8. 1923058E+01 » 8 „ 88 1 508 7 E+0 1 , 9 . 6499 13 3E+0 1 , 1.0509380E+02 » 

2 1.1473976E+02, 1 .2360462E+02 , i .37889282+023 ' 1 .3:18358 1E+02 , ' 

2 1.677377 5 E+0 2 , 1 . 8595332E+02 , 2 . 0692 2 74E+02 * 2. 3119103E+02 , 

2 " 2 *3'94383'8E + 02 » ' 2 .22521 05 E+02 , 3 . 3X52694E+02 »' 3 .77332 J3E+02 , 

2 4.3330 26E+02 , 5 . 00277 1 4E+02 , 5 . 8192238E+02 , 6 . 825 197 8E+02 * 
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2 8.07"3374E*+02i ~9.6'6~3~54l4E+02 . 1 . r£94491 E+03T~1 . 43422 07E+0 37 

2 1.7860852E+03, 2 . 2641 687E+03 » 2. 93Q7839E+03 , 3. 8892080E+03 » 

2 5.3189242E+03, 7. 5502829E+03 , 1 . 1236044E+04*' 1 • 7787769E+04 * 

2 3.0642697E+04, 5.9664177E+04* 1 . 4098798E+05 » 4. 7435606E+05 » 

2 3.7T830331E+06/ 

DATAUF AU+DT+31T 5 D/ 1C570 OOQE+Ol> 1 . 5714186E+0 1» 1 . 3 7 32762E+01 » 
2 1.5763777E401, "1.58-073 15E+01, I 588"34"93E40iri . 393746 3E+0 1 » ~ 

2 1 • 60 144 1 3 E+0 1 s 1 .6109566E+01 » 1 .6218 185E+Q 1 » 1 . 6340569E+01 » 

2 1 . 6 477061 E+0 1 •» 1 ,.6 628 0452+01 » 1 .6793950E+0 1 *” 1 . 697525 2E+G 1 * 

2 1.7172478E+H1, 1.7386207E+01 , 1 • 761707 6E+0 1 » 1 . 786578 1E+0 1 » 

2 1.8133084E+01* 1 „ 84198 16E+0 1 , 1 . 8726884E+0 1 » 1 . 90"5327'4"E+0 1 * 

2 1.9406060E+U1 > 1 • 978041 1E+0 1 , 2 . 01.79599E+0 1 * 2 . 06050Q7E+0 1 » 

2 2.10 58140E+U1, 2'.I540637E+01* 2 . 20*5428~0E+01V~2 . 2W1014E+01 . 

2 2.3182956E+U1 , 2 .3802416E+Q1 , 2.4461916E+01 , 2 • 5 1642 14E+0 1 > 

2 2. 5’9 12"323 E-H) 1 9 2 .6T09545E+01 » 2 . 7339497E+0 1 » 2. 8466152E+01 » 

2 2 . 9433873 E+Ul » 3 .0467462E+01 » 3. 1572213E+01 » 3. 2753965E+01 • 

2 3 • 40 1 97 7 5 E+L) 1 > 3 , 537499ZE+01 > 3 . 6829344E+0 1 * 3 .8 39 1041'E+Q 1 » 

2 4.0069889E+01 » 4. 1876827E+01 » 4.3824079E+01 » 4. 5925337E+01 » 

2 4.819'5974E+U1 » 5 . 06"?3"28 1 E+0 1 , 5 . 33T6765E+01 »' 3 . 6Z0848IE+0 1 , 

2 5.9353433E+01 , 6. 2780045E+O1 , 6 . 6520729E+0 1 » 7 • 0612547E+0 1 / 

DATA(F9A'{ D *1 = 60*100)/ 

2 7.5098017E+Q1 » 8 . 0026082 E+0 1 , 8 . 545327 1E+0 1 » 9 . 1445 1 19E+0 1 » 

2 9.807789 2'E : H:)1» 1 • 054407 l‘E+02 » 1 . 13-63815E+02*" 1 .22793502+02 > 

2 1 . 33 0 5 279 E+0 2 » 1.4458987E+02 * 1 . 576 128 1E+02 » 1 . 7237209E+02 * 

2 ' 1C8 "1710'EE402T 2.0837'9’56E+02» 2 . 3 04 5 1 4 62 + 0 2 V'2 . 5-574 81 22+02,' 

2 2.855693 1 E+02 » 3 . 201949 1E + 02* 3 . 6094160E+0 2 ♦ 4. 0924100E+02 , 

2 4.66 48732+02, 5 . 3649902E+02 , 6 • 2X12740E+02 i 7. 2519715E+02 » 

2 8.54687 44 E + 0 2 , 1 .01793 2E + 03, 1 . 2268227E+03 * 1.4986157E+03 ♦ 

2 1.859 1 183 E+0 3 » 2 . 3480042E+Q3 , 3 .028 390 5E+03 * 4. 0047739E+03 » 

2 5.4585626E+03* 7 . 72329 1 3E+03 , 1 . 1457308E+04. 1 . 8082732E+04 » 

2 3. 105'880'3'E+04 » "6.0 Z01284E+04, 1 .420981 8E+05 i' 4.768052TE+05 , 

2 3.7926419E+06/ 

10 FORMAT ( 9 ( El 4o 7 ) » / ) ~ * 

11 FORMAT (415) 

12 FORMAT (8F10.+) 

13 FO-MAT ( 54H THE AMPL+TUDE OF -ES-ONSE = , 

1 E14.7 , / ) 

14 FORMAT ( 54H PHASE ANGLE DIFFtRENCE (DEGREE) = , 

1 E 1 4 • 7 »/'//) 

15 FORMAT l 1H 1 ) 

16 FORMAT ( 54H THE GUESS-ED AMPLITUDE OF' RESPONSE « » 

1 E14.7,/) 

17 FORMAT ‘ ( 5 ( E 14 . 7 » 6 X ") » /') 

18 FORMAT ( 5X6HQ ( LB ) = » l’/X3HHO = ♦ 12X8HM ( SLUG ) = » 14X6HD ( IN ) = » 8X1 2HW ( RAD/ 
1 SEC) = ) 

19 FORMAT (F10.5»7(E14.7»1X)»/) 

20 FORMAT T2X3TTHA = »~5K »'6HF I ( n*»9X.6HFl ("3 ) = » 9X3HFA4 , 12X» 6HFT ( 7 ) = ,9X » 3H 
1DA=»13X» 4HDA0 = »11X»2HF= ) 

21 FORMAT' T54H THE “C1TAR'ACT2—I5T IC FRCQUENCT » 

1 E14 • 4 » / ) 

22 FORMAT t54TT THE NONDTMENS 1'ON'AL I ZED FREQUENCY » 

1 E14C4,/D 

23 FO-MAT ( +5 » 5~X , ( 7FT0 • 5 ) ) 

24 FORMAT (54H THE ZERO-ORDER AMPLITUDE OF RESPONSE , 
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‘ "l"El4i7Y/)' " ' ™ 

25 FORMAT (54H THE NOND+MENSIONALIZED UPPER RESPQNSE= * 

1 E14.7,/} 

26 FORMAT (54H THE NOND IMENS 10NAL I ZED LOWER RESPONSE= > 

1 El 4 . 7 j /-') * ~ 

BKS=C1» BCS=C2, AN1=N1=2 • 5 » AN2=N2=2.5, TOL IS THE TOLERANCE OF 
ERROR FALLOWED FOR THE SOLUTION * BMAS5 TS " I'N SQJG, “AND D'E LT A TS 
STEP IN INCHES. 

AG ( I ) » AND AOC-U) ARE - I SERIES OF THE GUESSED VALUES OF AMPL I TUDES » 

A AND AO RESPECTIVELY. QN ( I ) THE GUESSED VALUES OF THE DYNAMIC 
LOAD IN LBF . WN"( I ) ARE SERIES OF THE NONDIMENSTONAL I ZED FREQUENCY 
TO BE USED IN THE CALCULATION. 

” -E'AD 12 » BKS'j BC'S » ANT," ATTZ ' ~ ~ 

-EAD 12, HO, bMASST DElTAT TOL 
READ 23, IG, I AG'( I) » I *1 * IG ) 

-EAD 23, IG, ( AOG ( l ) » I =1 »+G ) 

READ 23, 'LQ, (0NIT),I = 1,LQ) 

I A=0 . 

DO' 110" Tq=1',l0 *“ " 

-EAD 23, LW, <WN( I ) ,I=1,LW) 

P'RINT 15 
•f AA = 0C 
gS=QN(IQ) 

AMASS=BMASS/12. 

' ARJ DELSD'EL TA4WASS 
H0N1=H0**AN1 
H0N2 = H0**MAN2 
HOPNl=HONl*HO 
" HOPN2=HON'2*'HO 

WS2=BKS*AN1/AMDEL/H0PN1 

WS“SQRT( WS2 i ' 

P-+NT 21* WS 
DO 1TO +V/31 , LW 
+A=+Q~0 . 

wb=wn( rw> 

W-WB*WS 

PRINT 18 •• “ " “ 

PRINT 17, QS, HO, BMASS , DEuTA, W 

PRINT 22, "WB 

WB2=WB*WB 

WO=WS2*AMDEU 

Q=QS/W0 

" H O'M A X= H "0 * 0 '« "9~9 ’ 

PI=3. 1416 
BK=BXS/W0 

8C=BCS*WB/WS /AMASS 
PHO 1 = PT *H'0"N 1 
PH02=PI*H0N2 

BI=B"K7PH0T * ■ 

B3=BC/PH02 

THE A'MPi-ITUD'F OF A FOR W=0.» A0=0. IS ESTIMATED BY* LINEAR I NTERP05TT ION 

IF (W.GT.O. ) GO TO 50 

IAA=1 

AWO=Q*PI*HONi/BK 



* ~ TF '( AWO-F'l AT 5'G “43 

42 A=0.5*H0 
A0 = 0. 

GO TO 50 

‘41 n-1 ----- 

+2=50 

— CO" TO' ‘44 

43 +1=50 

1 2'= 10 0" 

44 I 3 = t 1 1+121/2 

" ' rD=i2~rT+ovr 

+F ( +D-1D 45, 45, 46 
4TS IF ‘ ( AWO'-F 1A" (T2T ) “ 47 49 , "48 

47 1 2= I 3 

“ GO' TO ' 44 ‘ 

48 11=13 

' " GO - TO" 44 ‘ ' ' 

45 AI=+1 

" - A=TAI + FAWO-FlAn ID ) /'f F1*ATI 27-FlA(+T)T) /iw.* 
A0=0. 

GO- “TO 50“ 

4 A=+3 / 100 , MHO 

A0=0. “ 

50 CONTINUE 

**“ GTE S'S' “TH E "V AID FS‘ “OF‘~A~~AN D“ A'O . ' 

IF (IAA-1) 1,6,6 
"T CONTTNCJE' ’ 

+A=+A+1 

f F - ('Ia.'GTTGT GO TO T1 0 

A=AG(IA) 

^ 0 -=- a - 05Ti - A , 

6 CONTINUE 
" — +N T 1'6T A 

I TER=0 • 

103 CONTTNOE' 

HA05={ l.+AO/HO)**2.5_ 

+A"07= ( 17+A07'HOT^F3% 5*H0‘ 

A+03A/U+0+A0D 

hasaho'+too; - 

HA=HA+1 * ‘ __ 

IHA=FA"‘“ ' 

4* A 4- 1 

(FI A^C THAI r-TTATIRK) V 

F+ ( 2D3F2AU+HA } + { HA-+HA ) #UF2AU++A1 ) -F2A ( I HAD } 
■ F 1 C3 r=T3‘A“C IFTA r + HTA- IFi'A ') * < F3 A U H AT ) -ft A( I HAD ) 
F I ( 4) =F4A { +HA ) +U+A-IHA ) * ( F4A ( IHA1 ) ~F4A( IHA ) ) 

F+1 5T3F'5'A'r-PH‘A‘)'+ {TTA-+H'A') *UF5 A ( IHA ID-FFAT I HAD ) 

FI ( 6 ) =F6A ( IHA } + ( H A- 1 HA ) * ( F6 A { I HA 1 > -F6A { I HA } ) 

FT'T7 T= FT A ‘TTH‘A > +T FTA-TR A )' F CFTA U H A 1 ) -F7A7 THAT)' 

FI { 8 )=F8A( I HA ) + ( HA- 1 HA } ■*• ( F8 A ( IHA1 J-F8AC IHA ) ) 
* — FID 9T=F9“A"{ THAT+7 HA-TFTA FF ( F9A ( IHA1)-F9‘A( IHA) ) 
F I ( 1 ) = FI(1)/HA05 

FT T2T= “FDTZ )' /H AO 7" 

F I ( 3 ) = FI ( 3 ) / HA05 
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- F IT 4 ) = TT( 4 ) 7HA07" 

FI (5J--FI (5)/HA07 
FT ( 6) =-FI'( 6 ) /hA0"7 
FI l 7 ) = FI <7)/HA05 
F r c 8 ) = FI('8)/HA07 
FI ( 9 ) = — FI < 9 ) /HA07 

F= ( FI ( 1 }*‘B1-A"*WB2T*'*"2h : ( F+ ( 3 ) *B3*A ) **2-Q*Q 
FA1=(FI I 1 }*B1-A*WB2)*{FI ( 2 ) *6 1-WB2 ) *2 . 

FA2=(FI'<3 )*B"3*A)*CFI C3)*B3+FI (4)*B3*A)*2. 

FA=FA1+FA2 
G = FI ( 7 ) — 2 • *P’I 

FA0=2.*(B1*FI ( 1)-A*WB2)*B1*FI < 5 ) + ( B3*A ) **2*F I ( 3 ) *F I ( 6 ) *2 . 

DEL®FA*n(9)-FAO*FT(8> 

0A=(-F*F+( )+G*FA0)/DEL 

DA0= ( F*F I (8')-G*FA )/DEL 

IF (ABS(DA) .G1 .0.5) GO TO 1 

A0= AO+DAO 

A=A+DA 

P-+NT 20- 

PRINT 19, HA, Fill), F I ( 3 ) » FA, F I { 7 ) » DA, DAO, F 

PRINT 10, FAQ, DEL « G 

ITER=ITER+1 

IF (A.L-E-.O.) CO TO 1 

AOB=ABS( AO) 

AAO=A‘B~S FA-A'OB ) 

IF t AAO.GT.HOMAX) GO TO 1 

IF { ITER. GT. 13) GO TO 1 

IF ( ABS(DA) .G1 .TOL) GO TO 103 

IF (ABS('DAO) .GT.TOL) GO TO 103 

R=F I C3)*BC*A/PH02/(FI ( 1)*3K/PH01-A*WB2) 

’ ALPHA=ATAN t R ) 

ALPHA=At-PHA*160. / 3. 1416 
P-RINT 24, AO 
PRINT 13, A 
AU=A+AO 
ADOWN=A-AO 
PRINT 25‘, ATI 
PRINT 26, ADOWN 
PRINT 14, ALPHA 
+AA=+AA+1 
10 CONTINUE 
END 
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Card 1 Format (8F10.5) 

BKS - Value of c- for the stiffness of the gas film force 

in lb £ . 

BCS - Value of c„ for the damping of the gas film force in 

lb f /ips. 

Ml - Power n^ for the stiffness force in terms of gas film 

thickness . 

M2 - Power ^ for the damping force in terms of gas film 

thickness. 

Card 2 Format (8F10.5) 

HO ♦ Normalized gas film thickness at equilibrium. 

BMASS •• Mass of the step ring in response in slug. 

DELTA •• Step depth of the pad in inches. 

IOL - Convergence factor for the amplitude iteration. 

Card 3 Format (15, 5X, (7F10.5)) 

IG - Total number of the guessed amplitude A of the 

response. 

AG(1) - Array of the guessed amplitude A of the response. 

Card 4 Format (15, 5X, (7F10.5)) 

IG - Total number of the guessed amplitude Ao of the response. 

ADG(I) - Array of the guessed amplitudes Ao of the response. 

Card 5 Format (15, 5X, (7F10.5) 

LQ - Total number of the force excitation to be investigated 

in the production run. 

QN(I) - Array of amplitude of force excitation in lb^. 

Card 6 Format (15, 5X, (7F10.5)) 

LW - Total number of the normalized forcing frequencies to 

be investigated. 


HN(I) 


Array of normalized forcing frequencies 
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If the sixth input statement is located within the do loop of 
IQ = 1, LQ, IQ sets of card 6 are required. 
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PROGRAM RSRK 1 T { I NPUT ♦OUTPUT * PUNCH » T APE99 ) 

COMMON Y » DY » ATABL, RTABL » IFVD* X, DX » W »" BK. » "T3C » Q, AMASS* DELTA 
COMMON HO* TLAST » AY* FAC* WB2 , HON1* HON2* HOPN1* HOPN2 » AN 1 * AN2 
COMMON XA> TA > XTA » IA, IA2, DX I ♦ WLB * Wd, oLOPE, TLW » WBO » XU 
DIMENSION Y ( 2 ) » DY(2)» A T ABL ( 2 ) » -T ABL { 2 ) * WO-K(18) 

DIMENSION XA(IOOO)* TA(lOOQ), XTAUGOO) 

EXTERNAL. DERIV* CNTRL* GRAPH 

11 FORMAT (8F10.5) 

12 FORMAT (54H MASS OF THE StAL (SLUG) » 

1 E'l 4 » 3 » / ) 

13 FORMAT (54H DEEPTH OF STEP DISCONTINUITY (IN) » 

1 E14.3,/) 

14 FORMAT ( 54H NCND IMENSI ONAlIZED EQUILIBRIUM POSITION, HO » 

1 E 14* 3 , /) 

15 FORMAT ( 54H FREQUENCY OF PERIODIC FORCE, W (RAD/SEC) , 

1 E 14 • 3 , / ) 

16 FORMAT ( 54H AMPLITUDE OF THE PERIODIC FORCE APPLIED (LBS) , 

" 1 El 4 • 3 » /’) 

IT FORMAT ( 54H NOND I MENS IQNAlIZlD DAMPING COtFF ICIENT OF uAS FILM C , 

1 E14* 3 , /’) 

18 FORMAT (6UH NOND I MENS I ONAl I ZED MULTIPLE OF STIFFNESS DISPLACEME 
1NT, K»E14.3*/) 

19 FORMAT ( 54H NOND I MENS I ONAl 1 ZED DISPLACEMENT , 

1 E14.3’,/) ~ ■ 

20 FORMAT (54H THE MULTIPLE OF TIME AND FREQUENCY (RAD) » 

1 E14.3,A) 

21 FORMAT <3X7HDfcGREE=,4X, 1 1HPHASE ( RAD > = , 7X 13HDI SPL ACEMENT = , 11X9HVELO 
1 C ITY = * / ) 

31 FORMAT ( 7X ♦ 23HDAT A ( XA { I A ) » I A= 1 » 70 ) / *2 ( E14. 6 * 1H » ) ) 

'32 FORMAT ( 5X » 1K2 »"1X *4 ( ET4« 6 »1H * ) ) 

33 FORMAT (7X»22HDATA ( TA ( I A ) » I A= 1 , 94 ) / , 4 ( F9. 3 » 1H, } ) 

3'4 FORMAT ( 5X » 1H3 * IX » 6 (T9 . 5 * 1H'» ) ) 

35 FORMAT ( 7X26HDATA ( XTA( I AP ) * IAP = 1 ,360 ) / *2( L14.6, 1H » ) ) 

36 FORMAT T5XTH2 ,TX'»4 { El"4.~6 » 1H , ) ) 

40 FORMAT (52H ERROR RETURN* DX = 0 ,/) 

41 'F'ORM'AT (“5 2H'~NORM'AL RETURN ~ ,/) 

42 FORMAT (52H ERROR RETURN* VARIABLE INTERVAL MODE ONLY »/) 

TLAST IS THE FINAL TTME OF THE INTERVAL TN INTEGRATION* 

TLW IS THE FINAL TIME THAT THE SMOOTH VARIATION OF FREQUENCY ENDS. 

WLB IS THE TTNAL NOND IMENS TONAL I ZED FREQUENCY TO BE USED IN THE 
TECHNIC OF SMOOTH VARYING FREQUENCY. 

Yd), -YT2), X' AR'E INITIAL CONDITIONS OF X* DXTDT ♦ T RESPECT TVE'LY . ' 
BKS=C1. BCS=C2. DELTA +S STEP IN INCHES. QS IS DYNAMIC LOAD IN 
LB. A'M TS MASS 'OF THE PAD IN SLUG." WB IS THE NOND IMENSI ONAL IZED 
FREQUENCY. AN1=N1=2.5, AN2=N2=2.5 DXI Io THE TIME INCREMENT IN DEGREE. 
READ IT* TLAST, TLW* WLB 
READ 11* Y( 1) » Y ( 2 ) » X 

READ 11 * ETfCS * BC‘S, 'QS, AM, DElT'A, HO, W8‘ 

0+ 

H0N1=H0**AT11 

HON2=HO**AN2 

HOPNl=HONT*HO 

H0PN2=H0N2*H0 
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AMAS$=AM/121 

AMDEL=AMASS#DElTA 

WS2=BKS*AN1/AMDEL/H0PN1 

WS=SQRT(WS2) 

W=WB*WS " ’ ' 

WB2=WB*WB 

WL=WLB*WS 

XO=X 

SLOPE=( WLB-WB5 /WB/( TLW-XO) 

WB03WB 

WO=WS2*AMDEL 

Q=QS/WO 

BX=B'KS/WO * - - - - 

BOBCS/WS/AMAES 
PRINT 12, AM 
PRINT 13, DELTA 
PRINT 14, HO 
PRINT 15 , W 
PRINT 16’; STS ‘ 

PRINT 17, BC 
PRINT 18, BK 

DEFINE THE INITIAL VALUE OF Yd) 

XAU)=Yn> 

TA { 1 ) =x 

XT A ( 1 ) = Y T2 ) * 

DEFINE X AS Yd), AND DX/DT AS Y [ 2 ) 

THE FOLLOWING IS PREPARED FOR RUNGE— tCUTTA NUMERICAL ‘INTEGRATION* 
PI =3* 1415926535 
DX = DXI#PI /r3o* 

FAC= 10. **0 • 2 

NTR'Y=i ■ ..... 

N=2 
I FVD= 1 
I BKP= 1 

AT ABL ( 1 ) -0 • 0 01 
AT ABL ( 2 ) =0 • 0 01 

RT'ABL ( 1 ) =0*00] ' ... 

RT ABL ( 2 ) =0 • 001 
IA2 = I 
I A= 1 

‘ PRINT 21 ' ‘ 

CALL RKS3(DERIV,CNTRL»Y»DY»ATAbL,RT ABL , WORK. »X»DX*N»IFVD»I BKP »NTH Y 

1 IER'R) 

IF (IERR-u) 4, 5, 6 

4 CONTINUE' ' ... 

PRINT 40 

GO T‘0 3 " " " 

5 CONTINUE 

' PRINT 41' .... -- - 

GO TO 3 

'6 CONTINUE " ... 

PRINT 42 

3 CONTINUE " - 

THE FOLLOWING IS FOR THE PLOTS OF X AND DX/DT V.S. T. 
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CALL NA'MPlT 

CALL SCALE l XA , 5 . 0 » I A , 1 ) 

CALL SCALE { TA , AO . 0 , I A i 1 ) 

CALL SCALE < XTA » 5 . 0 , I A , 1 ) 

CALL AXTS(0.0,0.0,7HT VALUE i 7 * 40.0 » OTO »T"A ( 1 A+lTiTAT T‘A+2 ) ) 

CALL AXIS ( 0 *0 » 0 , 0 » 7HX VALUE *-7 » 8 . » 90 , 0 » XA ( I A+ 1 ) ,XAC I A+2 ) ] 

CALL AXIS { iVC » 0 • 0’,8HXT VALUE, -8 ,8.0790.0 ,XTA(TA+I FiXTA ( TA+2 ) T 
CALL LINE (TA,XA,IA. 1,1,0) 

CALL LfNE t FA » X T A , I A , 1 , 1 , 3 ) 

CALL SYMBOL ( 3 .0 ,9#0 ,0. 20 » 16HPL0T OF X V.S. T, 0.0, 16) 

CALL SYMBOLIC, 0,10. 0,0. 20,2'OHPlOT OF ‘DX/OT V.S". TfOiO»15J 
CALL ENDPlT 

THE FOLt-OWFNG IS FOR THE PHASE PLOT OF" X V.S. DX/DT.' 

Rl A= I A 
WIDTH=1Q. 

WSPACE=10. 

CALL S’ETUPC ( 0 • ,0 . ,0 . , 1 . , -20. , 2 . , W I DTH , WSPACE ) 

CALL PARAM ( GRAPH » 1. ,1. ,RIA) 

CALL ERDSURF 
END 


OB-OUTTNE DERI V" *' ' 

DIMENSION Y(2), DY t 2 ) * ATABL ( 2 ) , -TABL12), W0RKC18) 

DIMENSION XA { 1000 } , TAtlOOO), XTA(lO'O'O) 

COMMON Y, DY, ATABL, RTABl, IFVD, X, DX,W, BK, BC, Q, AMASS, DELTA 
COMMON HO, TLAST, AY, FAC, WB2 , H0N1, H0N2, H0PN1, H0PN2 , AN1 » AN2 
COMMON XA, TA, XTA, IA, IA2, DXI, WLB » WB, SLOPE, TLW , WBO, XO 
DYX 1)=DY( l’l /DT, OY'<2T=DY(2)/DT. 

DY(1)=Y(2) 

DY ( 2 ) = ( BK./H0N1 +Q*COS ( X ) -BK./ C HO-Y Cl)) **AN 1-BC^YT 2 ) *WB/ C HO-Y CD) **AN 
12 ) /WB2 
RETURN 
END 


SOBROUTINE CNTRL(NTRY) 

DIMENSION Y ( 2 ) » DY ( 2 ) , ATABLC2), RTABLC2) 

DTMEN'S ION XAT1000), TA~( 100U) , XTAC1O00) 

COMMON Y, DY, ATABL, RTABl, IFVD, X, DX»W, BK, BC, Q, AMASS, DElTA 
COMMON HO, TLAST, AY, FAC, WB'2 » H0N1, H0N2, H0PN1, H0PN2,AN1 »AN2 
COMMON XA, TA, XTA, IA, IA2, DXI, WLB, WB, oLuPE, TLW, WBO, Xu 
51 FORMAT' ( 5X » I 5,4( E"14. 6 »6X ) ) 

54 FORMAT (54H THE DISPLACEMENT IS OUT OF RANGc. ) 

IF YCl) IS' LESS THAN “4. » OR Y(l) IS~ GREATER IB AN HO, TERMINATE 
THE NUMERICAL INTEGRATION. 

A’Y=ABS (Y(l) ) 

IF (AY.GT.4.0) GO TO 2 
TF (YCl) -HO)' L, 2, 2 
2 CONTINUE 
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PRTNT 5'4 
NTRY=2 
1 CONTINUE 

IF (X-Tl-AST) 3, 3* 4 
4 NTRY-2 
3 CONTINUE 
I A2='I A2+1 
I ADX=IA2*DXI 
I A = I A 2 /4-F 0 »'0‘0 0‘1 
XT A ( I A ) =Y ( 2 ) 

XA(IA)=Y(1) 

TA ( I A ) =X 

PRTNT '51 '» I A'DX » X, YT T) > Y*( 2 } » WB 
IF (WB-WLB) 11, 12 » 12 

11 WB=WB'0*(1.+SL0PE*(X-X0') ) 

GO TO 13 

12 WB=WL6 

13 CONTINUE 
WB‘2 = WB#WB 
RETURN 
END 


SUBROUTINE GRAPH l T *XP » YP >ZP ) 

COMMON Y1 DY > ATABL » RTABL » IFVD» X»~'DX»W* BK. »”BC » "Q * AMA'SS» DELTA 
COMMON HO > TLAST, AY* FAC» WB2» HON1* HON2* HOPN1* HOPN2 *AN1 » AN2 
COMMON XA» TA s XTA» TA* IA2* DXI* WLB » WB» SLOPE* TLW * WBO» XO 
DIMENSION Y ( 2 ) » DY { 2 ) > ATABL(2)» RTABL(2)» WORKI18) 

D I MEN'S I ON XA(IOOO), rAllOOO)* XTA(IOOO) 

NP=T+.0001 

XP=XAfNP)*5. 

YP=O.C 

ZP=XTA(NP)*5. 

RETURN 

END 
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Card 1 Format (8F10.5) 

TLAST - Final normalized time in radian for the integration. 

TLW - Final normalized time in radian for the changing of 

Tf during the integration. 

WLB - Final normalized forcing frequency. 

Card 2 Format (8F10.5) 

Y(l) - The initial value of X. 

?(2) - The initial value of X, i.e. dX/dT. 

X - The initial time, T^ for the integration. 

Card 3 Format (8F10.5) 


BKS 

- 

Value of c^ in lb^. 

BCS 

- 

Value of Cj in lb^/ips. 

QS 

* 

Value of q in lb^. 

AM 

- 

Mass of the ring in response in slug. 

DELTA 


Step depth, 6 in inches. 

HO 

- 

Normalized equilibrium gas film thickness. 

WB 

- 

The initial normalized excitational frequency 


Card 4 Format (8F10.f) 

AN1 ■- Value of n^ which is 2.5 in the case being investigated 

here. 

AN2 - Value of n„ Which is 2.5 in the case being investigated 

here. 

DXI - Time increments during the fixed interval of the numerical 

integration. 



PR0GRAM“T7STAB~C INPUT » OUTPUT ) " "" 

DIMENSION SA ( 50 ) » SB(50), SC(50)» SD(50)» SE l 5G ) » SFC5Q) 

DIMENSI ON PO( 50) >' PCRESTH » PC I ( 50) » "OX l 50 ) » C(TOO)* *170*70) 
DIMENSION CI(2»3»50)» SUBVI50), VI(20) 

' DIMENST ON" XXCIO'O ) ♦ 8 C5UV3 ) ♦ BH ( 50 ] » AH ( 5 0 ) '» CHI 50 ) 

COMMON B I » AEAM * HI* T0L» ML » N2* B2» NL 1 » N2l» NP* NP1* NL2* Nl3 

” ~ COMMON" TI2V H 13 i 'H'23 » "HSUM ' “ " — 

COMMON SA* SB* SC* SD> SE» SF 
“li form'At" ($FrorA;3T5) ~ ~ ■ ■ — ■ 

12 FORMAT (8F10.5) 

13 FORMAT”! 5'4H'"THE THRESHOLD “FREQUENCY” OF THE STEP~5E'ATr * 

1 E14.4,/) 

~*14 FORMAT T5A FTTH E S'Q'UE EXE' “FTl'M PARAMETER * 

1 E14.4,/) 

'T5 FORMAT T54H THE BEARING "NUMBER ' » 

1 E14.4,/) 

1"6 F'ORMAT" (54H MASS " » 

1 E14.4,///) 

’IT FORMAT T54TTRATI0 “OF" MITT CLEARANCE HETGHT'TO DTFF.' TN CL'EA'RANCES » 


1 E14.4,/) 

IE FORMAT C54H THE 1ST ORDER PRESSURISED' LOAD -- - t 

1 E14.4,/) 

19 FORMAT (T8Et4-.4) ,/-) ~ 

20 FORMAT ( 54H THE POINTS ARE TOO FEW ) 

”21 FORMAT" (SAH'A' AND/OR EFTS 'OUT 'OF RANGE -OFT 'ABLE" * “ " ) 

22 FORMAT (54H THE GUESSED THRESHOLD FREQUENCY OF THE STEP SEAL » 

*“ 1 El"414» D - - - 

23 FORMAT (52H THE GRID DIFFERENCE ,/) 

~”24 F0RMAT T32H"TKE rERO~0RTYER PRESSURE"'DTSTRTtfUTTUN ” ,/) 

25 FORMAT (54H THE VALUE OF SMALLNESS , 

" ' 1 EI4V4T/) * — 


26 FORMAT (54H LOAD DUE TO THE IMAGINARY PART OF COMPLEX PRESSURE * 

• 1 ET 4 • 4V7") • 

27 FORMAT ( 5X * I 5* { 7F 10.4) ) 

“"29 FORMAT" K7 *~C8‘E14«'4'»7 ) ) '' ” "" ----- 

30 FORMAT I52H THE IMAGINARY PART OF COMPLEX PRESSURE DISTRIBUTION*/) 
3"1 FORMAT ( 5"2H "THE' "REAL "PART OF COMPLEX" PRESSURE'DT STRTBUT ION ' ,/) 

32 FORMAT (54H THE RATIO OF WIDTH OF STEP WITH HEIGHT HI , 

I ET4V4 * / ) " ------ - ' 

READ 11* Bl* AlAM, HI* TOl, NL, N2* NV 

" RE A'D" ' T2 > TVnTr*T = l*'N\T) “ " 

B2=1.-B1 

N'Er=NCTE — — -* 

NL2=NL-2 

" ' NL3=~NL-3'--' — • ~ ' 

ND=NL2+NL2 

N2T=NX-r - ... 

NP=N2+1 

-npi=tip+t — ’ - ---- — 

H12-H1**2 

H13‘=T1W3”- ' 

H2=H 1+ 1 • 

"H2T=H2*'*2 ' ' ■ ■""* 

H23=H2**3 



' "READ ' 1'2 V TB’XTTFiT =1 ,’NUiT ' “ “ * ~ 

READ 12, (POII) ,1=1, NL) 

RSQW-Hl3T+3«’ : S r Hl¥H2'^‘(Hl-< : H2 )+H23 

CALL NEWR<PO,DX) 

CALL" TN’TEG’t O'* » IV »'DX »PCT , NL'» WLOAD , C I , TERR ) "" 

PRINT 32 , B 1 

“'PRTNT~TT, ~AUAM' — 

PRINT 17, HI 

PRTNT-pr, ~srG~ • ‘ ■ ” ■ 

PRINT 23 

prtint T9 ,'■■( Dxrrr,'T= r^nr ■ - - — — - 

PRINT 2 A 

* -pRTNTT^T“(wm“r=-"iTmrr — *'~* ~~ — r — 

PRINT 25, TOL 

' _ Al“=-ALA7T**HT72‘ a - 

A2=— AlAM# ( H1+H2 ) / A# 

A3 h=Al"A'M¥K2T2V ■“ ■* - ” 

THE COEFFICIENTS OF SMALL A* B, C, D, E» F AT EACH GRID POINT ARE 

'CATCITCATED' FN T HE~F OUCOWTN G~'A'S 5A, .T. ETC*' ”* 

DO 2 1=2, N21 

SATIT=HT3/DX ( I ) ’ ~ ' ' 

set I ) =H13/DX { 3 — 1 ) 

■ sbti r=~sA 1 1 r-s'c rn • - — — 

SDl=(POt 1+1 ) **2/DX t I ) -PO( I )**2*( l./DXU )+l*/DX( 1-1) ) +PQ t I-l)**2 

1 TD'X'CI-n r¥3T^Tir27T. - ■ 

SD2=ALAM*(POt 1+1 )-POt I— 1 ) )/2. 

SDTD =S'D1-SD'2' ‘ “ ' ' 

SEl I ) =H1* t DX ( I)+DX( I— 1 ) j/2. 

SFTIT=TDX ( IT+DXn-TT'r/T. ' 7 ~ 

AH ( I ) = A 1 

BH( I ) =0 • 

2 CONTINUE “ ' 

SA(N2) =HSUM/8./DX(N2) 

SC t'N2T=F)137DX(7T2T )’ "" * ” * - 

SB ( N2_)_=-SA ( N2 ) -SC ( N2 ) 

SD1= t t POTNPT **2-P0 (772 )'** 2' )* t H 1 2 + ff 1 *HZ+ H2TT7D XT7T2T+ TPO CN ‘271** 2 
1 PO(N2)**2)*3.*H12/DX(N21> )/A. 

S'D'2= ALA'M* ( P‘0 (W ) -PO CN2T') ) 72 . " 

SD ( N2 ) =SD 1-SD2 

SET N2 r= R H T+H 1 ) TDX l N2T + 2 . *H V* 0 X ( N 2 1 Y T 74 • ~ 

SF(N2)=(DX(N25+DX(N21> )/2. 

AFTTnT )"= A2 “ ' ~ ■ ■ 

CH (N 2)=-A1 

BH ( N2 ) = AH ( N 2 ) +CH ( N 2 ) ‘ “ ' ‘ ~~ * **~ 

DO 3 I=NP1,NL1 

~S ATT r='FT2'37 0XT T) ~ 

SCt I ) =H23/DX{ [-1) 

■ smTr==SATrr-scm — - — - — 

SDl=(POt I + l)*fr2/DX( I )-P0{ I )**2*< l./DX( I )+l*/DX( 1-1) )+P0( I-l)**-2 

■ 1 TDXn-Tn ~*T <f f H 227 4V ' 

SD2=ALAM*(PO ( I + D-P0C1-1) >/2. 

*SDriT=SDl=^SD2 •“ 

SE! I ) =H2* t DX t I )+DX( 1-1) )/2. 
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" ~ sttit=idx rrr-FDX'CT-n rrz', ~~ 

AH ( I )=A3 

■ cHCir=-A3‘ ■ ' ~ 

BH( I ) =0 • 

' 3 CUNT 1 NUF ' ‘ ' 

SA(NP)=H23/DX(NP ) 

S C TN P7 = H SO M7~8Y ! D XTN2 5 ' — ~ ~ "" 

SB ( NP ) =“SA ( NP ) -SCCNP) 

SD1* r ( FO fNPTT* *2 - P 0 { NPT**2) *3 .*H 22 /DX~( NP r+TPO l N2T**2=POTNP ) **2T* 
1 (H12+Hl*H2+H2.2) /DX(N2) ] /4. 

' ' SD2=ALAPI*'lP0'CNPiy-P0CM2T')'/2. *" 

SD(NP)=SD1-SD2 

' “ S'ETNFT=T'{TT1TR2 ) * DXT N ZT+ 2' . *H 2 * D XTNP ) T T4\ 

SF ( NP ) = ( DX ( NP ) f DX ( N2 ) ) /2. 

■” AHlNP } = A3~” - - - ~~ 

CH(NP)=~A2 

BH ( NPV^ARTNP )+CH'CNP ) ------- — 

CALCULATE ELEMENTS OF (A) AND (C) IN THE EQUATION AX=C. 

“■*" B NfATRTX'TS“~USED TO SAVE - AH *'J ) WH fCRTARE' TN DEFEND ENT ON "THE 
DO 4 J= 1 » ND 

~ DO" 4 I'=T»ND “* • 

A ( I * J ) =0 . 

'4 CONTINUE ‘ ” - 

V-V I l 1 ) 

' ■' P'R'INT ‘ 2T, ’V"'” “ ~ 

AMASS=WlOAD/ V*#2 

PRINT T6T "AMASS ~ “ 

DO 5 J = 2 » NL3 

I=J+I “ * ' ' 

JN=NL2+J 

'* ' TNT3NTT * "• — 

AC J* I }=SA( I )*POC 1 + 1 > +AH C I ) 

ATUN » IT5T=“A CJ * T ) " 

AC J»J)=SB( I )*POl I ) +BH C I ) 

ATJNVJNT=AT3'»J) ' ~ 

A(J*J-1)=SC( I)*PO( I — X ) +CH ( I ) 

ATO N“» UN ~T )-~ATJ TU^TT ' ----- - -- - - - — 

AC J*JN)=SE( I )*V 

A (ON >'JT=- A ( J”» JN ) “ " * ' 

C( J)=-SD( I ) 

CCCJNTXSFTIT-SFOC IT*V ’ “ 

B ( J » 3 ) =A( J » J+l ) 

BTUr2T=A~t JIUT' ' ” ‘ ' 

B £ J » I ) =A ( J » J-l ) 

5 CONTINUE - ----- . 

A ( 1 » 1 ) =S8 C 2 ) #P0 ( 2 ) +BH ( 2 ) 

- ATNcmrciT^AT i'm ■ - 

A c 1 9 2 ) =SA ( 2 ) *P0 C 3 ) + AH ( 2 ) 

~ "'ATNinTNTrS'All »'2T ‘ ~~ ~ 

A(1»NL1)=SE(2)*V 

“■ "A tn lt rn^'-ATi tnld - • — - - 

C(1)=-SD(2) 

C TN inT="SFT 2T * P QT 2') *• V“ ' * ' ' 

A(NL2,Ni-3)=SC(NL1)*POCNL2}+CHCNl1 ) 



A TND VttD'-l') = ATNV2 ;nm ~ * * 

A ( NL2 » Nt -2 )=SB(NL 1 )*P 0 (NL 1 ) +BH ( NL 1 ) 
A<ND»ND)=A(*Nl 2 ,NL 2 ) 

A(NL 2 »ND) =SE(Nl .1 >*V 
A ( Rt> > NU 2 T=-fi (NL 2 » ND* )' 

C < NL2 ) =-_SD ( NLl ) 

CTN*D)'=SF"(KLl )*“P'0'(¥L1 )*V~ 

8 ( 1 * 1 ) =A ( 1 * 1 ) 

" " bti * 2 r=m » 2'i 

B(NL 2 , 1 ) = A { NL 2 , NL 3 ) 

"" * B ( NL 2 » 2 ) = A ( NL. 2 , N L * 2 ) “ 

DO 61 J= 1 »ND 

xxtd ) =crjy 

61 CONTINUE 

CALC TSETE'QTA ,XX , ND'VD'ETT 

IF ( DET- 0 • ) 108 , 109 , 108 

“ro ‘9 STOP ' ■ ■“ 

108 CONTINUE 

p cT - (Tr= - 0 --- - - - " " 

PCI (NL 3 = 0 . 

pc r ( nw.* “ 

PCR(NL)= 0 . 

DO "7 “D = NLlTNO *~ * 

PC I ( J-Nt -3 ) =XX ( J ) 

~ ‘ 7 C'O'NTTNtTE “ ‘ " " * ~~ 

PRINT 30 

~~ P R TNT 2*97 * ( PCT (II » I = 1 TNL ) 

CALL INTEGt 0 . , 1 . ,DX,PCI ,NL,PSrtl ,C I » I ERR ) 
IF T Em - Y) 113," 114, "115 

114 PRINT 20 

STOP "~* “ " 

115 PRINT 21 
STOP 

113 CONTINUE 

PRINT 2 ” 6 ‘» “P 3 W 1 “■ ' 

DO 69 J = l, NL2 

— 6*9 PCR C3+rr=XTCJ ) * * * 

PRINT 29 -, (PCr-Un ,I = 1,ND 

P~R I NT “ 3 T - - - 

CALL I NT EG ( 0 • , 1 • » DX > PCR ,NL » PRSM » Cl » I ERR ) 
PRTNT'TBVTR'SM ‘ * 

AMASS=PRSM/V **2 
' PRINT - 1 * 6", "ATOSS ““ 

DO 121 I V = 2 »NV 
' v‘=VT ( fvT "*' 

PRINT 22 , V 
DO " " 5 T ' J=T»’MD*~ 

DO 51 1 = 1, ND 

ATr»”JT='0'; * ‘ 

51 CONTINUE 

DO 52 "J= 2 YNU 3 * ' 

I = J +1 

“JN-=NL2TJ-- 

IN=JN+1 
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ATJTJ- rr-BTTTT ) — — 

At J,J)=B( J,2 ) 

- * ' A( j*» J+ rr=BT( j,'3> - - — — 

A< JN*IN)=A< J»I ) 

ATU N * JNT = A ( J * J ") ' ~ ' — 

A( JN* JN“1 ) =A [ J » J-l ) 

‘"52 CONTINUE'” — * ~~ 

A ( 1 » 1 ) =B ( 1 * 1 ) 

‘ A (1 , 2 ) =BTI » 2T " • 

A(NL2»Nl-3 )=B (Nu2»1) 

A CNL 2 »'Mil2 f=BTN i_2 ,2 )' " * ' * " "" 

A(NL1»Nl.1)=AC1*1> 

'ATRLXT NU) =~AXT f 2 ] * 

At ND»ND-1 ) =A tNt_2,NL3 ) 

ArNT>;Ntn=ACNL2YNL2T ' “ 

CALCULATE THE ELEMENTS of A AND C WHICH VARY WITH RESPECT TO V. 
DO 6'J"-1,NL"2 ' 

I=J+1 

JfTTRNLT' “ 

At J»JN)=SE( I )*V 

“V“ AT J N~» 0 T= ~~A" < UTJNT 

C( JN)=SF( I )*PO( I )*V 

~6 CONTINUE- ------ - 

DO 62 J= 1 »ND 

- xxru ) ^ctot ~ - ~ - — " — - 

62 CONTINUE 

"CALL DETEQ(A'»X'X>ND»DET)' 

IF t DET-0. ) 111, 112, 111 

- 1X2 STOP' ------- 

111 CONTINUE 

~ ”* 'DO - ' 8 J='NLX,ND ' “ ~ 

PCI ( J-Nl.3 ) =XX ( J ) 

' 8 CONTINUE * 

PRINT 30 

PRINT" Z9T * t P"C1 ( rr»T=lTNL ) ~ " 

CALL INTEGIO. »1. »DX,PCI »NL*PSM2*CI, I ERR ) 

if 1 1 err- it" ir6nirrn8- - ~~ 

1 17 PRINT 20 

STOP ' — 

118 PRINT 21 

" STOP ‘ " ~ ~ 

116 CONTINUE 

PRTNT ”2’6V' PS>12 ~ ” — 

DO 9 J=l, Nl2 

9'"PC R (”J+ 1T=XX"CJ ) “ -- - - 

PRINT 31 

PRTNT ~Z9 V T PCR fT) '» I =T7M. ]' 

CALL INTEGtO. >1 . ,DX ,PCR,Nl,PRSM, C l , I ERR) 

“ P'RTN T ~ ITT* ""P'R'S M - 

AMASS=PRSM/V**Z 

‘ ‘ PRTNT “167 “ATTASS “ 

,121 CONTINUE 

FND - • ” ~ ‘ ~ 
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SUEfROUT TNE 'I NTEG'CAl 8YFT7F » NP , VALUE , C » TERR ) 


■ •*■)«•***#***■«-**-***■*****■*********•**-*■** INTEG002 


* I NTE60Q3 

SUBROUTINE INTEG *INTEG004 

• _™__ " - ‘ " *INTEG00'5 

♦INTEG006 

I NTBGRAT'ES THE NON 'EqDTDTSTA'NTlY T ABDICATED" 'FUNCTION' TTX ( I ) ) ” *INTEG00T 

BETWEEN THE LIMITS A AND B. *INTEG0Q8 

' *TNTEG009 

A MODIFIED METHOD OF OVERLAPPING PARABOLAS IS tMPLOYED. *INTEGQlO 

‘ ‘ * INTTEGOll 


A SECOND ENTRY POINT 1 LNTEG2* IS PROVIDED FOR MORE THAN ONE 
‘ INTEGRATION 'ON THE 'SATIE "DTV’T S IONS 'OF X. THT5 SAVES'THE TIME 
OF CALCULATING THE WEIGHTING FUNCTIONS. 


ARGUMENTS - 


"A 
B 
' 'X 


F 

NP 

VALUE 

C 

I ERR 


LOWER LIMIT' O'F ‘ INTEGRATION. ‘ 

UPPER LIMIT OF INTEGRATION. 

ARRAY OF ARGUMENT VALUES. MOST BE "MO NOT ON TCA L L Y 
INCREASING AND MUST BE DIMENSIONED NP. 

ARRAY OF “FUNCTION VALUES. MUST BE DIMENSIONED NP. 
NUMBER OF POINTS. NP MUST BE GREAT tR THAN 3. 
RESULTANT VALUE OF THE INTEGRATION'. 

WEIGHTING FUNCTION PASSED TO THE MAIN PROGRAM 
"FOR '"STORAGE. " ' ' - - 

RESULTANT ERROR PARAMETER. 


REQUIRED SUBPROGRAMS - NONE 


COMMON STORAGE - 

THE WEIGHTING' FUNCTION ~C I'S STORED I N* THE MAIN PROGRAM AND 
REQUIRES THE FOLLOWING DIMENSION STATEMENT WHERE D.GE.NP. 
DIMENSION C I 2" » 3 * D ) 


ERROR INDICATIONS - 
I ERR = 0 INDICATES NO ERROR. 

IER'R = 1‘ INDICATES' NP" IS LESS THAN 4. 

I ERR = 2 INDICATES THE LIMITS OF INTEGRATION ARE OUT OF 
THE RANGE OF 'THE TABLE. 
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*INTEG012 

*IWTEG0r3~ 

* I NT EGO 14 

* I NT EG 01 5 
*INTEG016 

* I NT E6017 
*INTEG018 

* ITTT E G 0T9 
*INTEG020 
*INTEG02T 
*INTcG022 
*INTEG02'3 
*TNT EG024 
*TNTEGQ25 
*T NT EG026 
*1NTEG027 
*INTE6028 
*1 NTEG029 
*INTEGQ30 
* I NTEG03T 
*INTEG032 
*INTEG033 ' 
*1 NT EG034 
*INTEG035 4 
*INTEG036 
*ITJTEG037 
*INTEG038 
* I NTEG039 
*INTEG040 
*1NTEG041 


*T NT EG042 


DIMENSION X ( 50 ) » F(50)» C(2,3,50) 


' DIMENSION H T 5 0' ) VSU BV ( 5 0") INTEG045 

INTEG046 

'* NP MOST BE GREATER' THAN' 3 ' INTEGG47 

INTEG048 

' IF FNPT EE'. 3 3 GO TO 96“ " “ INTEGOW 

INTEG050 

- •' CALCULA'TTON'“OF INTERVALS’ OF X INTEG05I 

INTEG052 

NH-NP'-T -■ ‘ I N'T EGOS 3 

X(D=0. 
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— D'crrcr t=itntt 

10 X ( I+l ) =x ( I ) +H { I ) 

DO '20' 1 = 1 »~N H ' 

IF (I.EQ.l) GO TO 15 

DEFINE COEFFICIENTS OF £IRSJ parabola^ __ 

C( 1 » 1 » I I) )**3/ 16.*H( I-1)*(H( I-l)+H( U n 
CT 1T2VD =R'( IT* T 3 .~*H a-T f+NCl ) ) / (6.*HTP-1 ) ) ‘ ' 
C(1»3»I)=H(I)*(3.*H(I-1) + 2 »*H ( I ) ) / (6**tH ( 1-1) +HK I ) ) ) 

15' CONTI’NUr 

IF (I.EQ.NH) GO TO 20 _ _ 

DEFINE COEFFICIENTS OJF SECOND PARABOLA 
C ( 2 » 1 » I ) =H ( I ) * ( 2 .*H{ I ) +3 ( I + l) )/ ( 6 •* ( H ( I ) +H { I + 1 ) ) ) 

cT 2 Y 2 »Tr=H"( n*(Tm )+3r*Hn+i n / (6**^ rr-f i ) r 

c ( 2»3 » I ) =-( H l I ) ) ** 3/(6. *H l I+l )* ( H( I ) +H (I + l))) 

20 CO'NTTNOE" ' ' 

— ’ ENTRY- INTEG2 


INITIALIZE SUMMATION VARIABLE 


VALUE=0 • 0 

" I'F' ' ( B-'A V 40 >92 » 30 

~ - 8 IS' GREATER’ THAN'TA 

"30 ALIM = 'A 
BLlM = B 
SI'GN =1.0 
GO TO 50 

A IS GREATER THAN B 

40 ALIM * B 

B'CIM = A 

SIGN =-1.0 

' 5"0 -NH = NP-T 


~S"ET TTNG" 7 H ETOV? ER “LTMIT 'OF rNTEGRATTO'N ' 

~B0‘ 63"T-T7NFT 
PARTA = 1.0 

IF { AL I M— X (T ) ) “ 6"1"» 6~9~» 6"3 
61 IF (I .EQ.l) GO TO 97 

ATTM ~— xTr-TH 

PARTA= ( X ( I J-ALIM) / tX( I ) -X ( 1-1) ) 

GO TO" '6*9“ ~ ~ 

63 CONTINUE 

— SR"' CONTINUE: ' ' 


1"NTE605“6 
INTEG057 
"INTEG0"58 
INTEG059 
“ ’ TNTEGO'6'O 
INTEG061 
TNT EGO 6 2 
INTEG063 
ITMTEG064 
INTEG065 
~~~ " INT EG 06' 6 
I NT EGQ67 
TNTEGO'68 
INTEG069 
TNTEG07Q' 
INTEG071 
TMTEG0T2' 
INTEG073 
LNTEG074 - 
INTEG075 
INTEG076' 
INTEG077 
TNTEG078 ' 
INTEG079 
TNT EGQ8Q 
INTEG081 
INTEG082' 
INTEG083 
' rNTEG-084 
INTEG085 
INTEG086 
INTEG087 
' - INTE6088 
INTEGQ89 
" ’ “TNT E'G09'0 
INTEG091 
I'NTEGO'92" 
INTEG093 
TNT £6094 
INTEG095 

- — ' "T'NTEG'0‘9'6" 

INTEG097 
T NT EGO 98' 
INTEG099 
" ~~ INTEGTOO' 
INTEG101 

TNTEGTOZ 

INTEG103 

TNTEGT04 

INTEG105 
' — IN TEG 106" 

1NTEG107 


£. £ t 
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SETTING THE TP'P ER~“L I*M1 T ' OF INTEGRATION ' TNT E'G’Fo'8" 


DO 73 I =1 »NH 
PARTB = 1.0 

IF' TBLI M-X (If) 71 .79173 ' 

71 IF ( I.EQ.l) GO TO 97 
' " BLlM = X( IT ” 

PARTB=(BLIM-X( I - 1 ) ) / C X ( I > —X ( I — 1 ) ) 

GO TO 79 “ • 

73 IF ( I.EQ.NP ) GO TO 97 


INTEG109 
1NTEG110 
INTEG111 
"I NT EG IT 2 
INTEG113 
TNTE6IT4 
INTEG115 
TNTEGT16 
INTEG117 


79 'C’O'NT I'NU E‘ ' ‘ ' * 

~ "‘C'AITCOC'AT TON 73F FRTE'GRTTL ' OVER SUBTNTERVAC ~ — 

DO ~8'0‘ 1 = 1 »NH 

SUBV( I ) =0.0 

I F T X ( I ) . EQ*. AL I M7 ' SUBV1T) =C( 2 » 1 . I ) #FT I ) +C 12 » 2VT)'*F n + 1 ) +C ( 2 » 3VI ) 
1FU+2) 

IF ' ( I -NIT)' ~TU 2 y T03V l'OT ” - • -■ - - 


TNTEGI18 
INTE6119 
TNT EG 12 O' 
INTEG121 
TNT EG 12 2 
INTEG123 
* TNTEG124 
INTEG125 


103 CONTINUE 


ADX = AB’S ( X ( I + IT-BLIM ) ' 


IF (ADX.LT.l.E-5 > SUBV ( I ) =C ( 1 » 1 » I ) *F ( 1-1 )+C( 1 »2* I ) *F ( I > +C( 1 .3* I ) * INTEG126 

1 FU + T) . - 7 NT E6127 ' 

GO TO 101 

T02 TONTITOJE" ‘ -™ 

IF ( X ( I) .GT.AlIM.AND.XII+1 ) .LT.BlIM) SUBV ( I ) =0.5* < CT 1 *1 * I ) *F ( 1-1 ) INTEG128 
■“ r+rcrr.2Yi )+~ct 2 »i» ny*TTi ) + <cu.3> i r+ci2»2"»i > t*ft i+r)+cr2“*3*i j* integtz9 - 
2FU + 2M INTEG130 

"101 CONTINUE" ' 


IF (PARTA.NE.1.0R.PARTB.NE.1) SUBV ( I ) =PARTA*PART B*SUBV U ) 

CALCULATE THE FINAL VALUE OF THE INTEGRAL 

80 VALUE=VALUE+SUBV ( 1) 

VALUEX5TGN* VALUE - - - 

SET ERROR PAR'AMET ER"“FTFr 'NORMAL RETURN ” 

92 IERR = 0 .... 

RETURN 

SET ERROR PARAMETER FOR TOO FEW POINTS 


INTEG131 
TNTEGT32 ' 
INTEG133 
I N T EG'l 34" 
INTEG135 
INTEGI36" 
1NTEG137 
‘TNT EGT37T 
INTEG139 
TFlTEGTAO"' 
INTEG141 
• I NT EG 142" 
INTEG143 
TNTEGT44' 


96 IERR = 1 INTEG145 

'RETURN " ‘ ‘ INTEGT41S" 

INTEG147 

~ SET 'ERR'OR'"P'ARA"METER _ FOR ' A 'AND/OR B OUT OF "RANGE" OF 'T AB1_ E TNTEGIW 

INTEG149 

”97 TERR"™ 2 ' TN'TE6T5~0' 

RETURN 1NTEG151 

” ' END - - - • I NTEGT52 
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SUBROUTINE D'ETEQ C A", B HT* DET ) ' " 

DIMENSION A (70, 70) »B(70) »IPV0T(70) 
DET = 1.'0* 

DO 11 J = 1 »N 
"T1 I P V OTTJ 1 = 0Y0 ~ 

DO 121 I = 1 » N 

~ ’■ T = 0.0' ” 

DO 51 J=1»N 

- -- x p “( TPVOTt J ) -1 P 21 > 5 1 *~2T ' ’ ' 

21 DO 50 K*1»N 

" " I F "I I PVOTC K ) -1) ' 3Y» 5 07 5~0' 

31 IF ( ABS ( T ) -A3S ( A { J » K. ) ) ) 41,50,50 

~7jrr TROW=Zf 

I COL=K. 

- - r=firovK') - 

50 CONTINUE 

~5T CONTINUE 

IF (ABS(T)-l.E-8) 131,131,55 

—55 ywotttcol ) =r " ~ ' 

IF ( I ROW— I COL ) 61,81,61 
•'"61 D'ET=-D'ET ‘ ~ ' “ 

DO 71 L = 1 » N 
T=At I R0 W7C) 

A ( IROW»u)=A( ICOL,L) 

“71 AIYCOCTlT) =T 

T = B ( IROW) 

B FI ROWT= B'( I'COL ) * ' 

B ( ICOL ) =T 

81“ T EMP-TCfl COT TTCO L T - * 

DET=DET*TEMP 
“ “ ATTCOL VI COLT=l . ~ 

DO 91 L=1»N 

91 A f ICOL ,1- 1= A T ICOL', L ) / TEMP" ' 

B ( ICOL )=B(ICOL)/TEMP 
DO’ 121 L1 = 1VN 
IF (Ll-ICOL) 101,121,101 
“101 T=A _ (ur, ICOL)* " ~ ’* “ 

A ( L 1 , I COL ) =0 • 

' DOTirX=I,N ' ----- 

111 A(L1»L)=A(L1,L)-A( ICOL,L)*T 
BTLT) =B'( L l )“-B (ICOL) *T“ 

121 CONTINUE 

' RETURN' ’ ~ “ 

131 DET=0.0 
"" RETURN' 
end 


SUBROO'TT N E“NEW RTP'OYDX) ' ' ' 

DIMENSION AA150), BB150), CC(50)» DD(50), E£(50)» FF ( 50 ) , DX(50) 
— TTMENSIUN' AT5 0~) » ~ BT50T, C( 50 ) , PO ( 5QT, DPT501 f F I C50 ) » FI DP ( SO - , 3 ) 
COMMON Bl> ALAM, HI, TOL» NL» N2, B2'» NLl , N21, NP, NP1, NL2, Nl3 



COMMON R 2 » ' HT3 ,"~TT2 3 , ~TTSU M ~ 1“ 

COMMON AA, BB» CC» DD* EE, FF 

11 FORMAT "(4F10«4,2I5) ' “ 

12 FORMAT (8F10.5) 

"TO FORMA T 1 6'X ETA* 4 ,'3 ( '6 X *ET 4.4'),/) 

14 FORMAT (17X3HFI«»3X17HDFI/DPIJ)»J«1»2*3*/> 

“T5 FORMAT ( 3X2HT=T22X3HDX^;T7X‘3HP0= V/ )' ~ 

16 FORMAT <5XI5,2(6XE14.4) ,/ ) 

‘17 FORMAT ("5~4H ‘PRESSURE Ts SMALLER T HAN'T HE AM b I ENT PRESSURE 

18 FORMAT (5XI5»3<6XE14.4) ,/ ) 

19 FORMAT 7 1HT) ‘ “ “ 

DO 2 1=2, N21 

" ” A~A7 n =H2372T7DX7T ) ~ ' ' 

CC( I ) =H23/2./0X( 1-1 J 

' ‘ “ BBT I ) =-A7\'< n-CC ( n ‘ " ' 

FF { I )=AuAM*H2/2. 

' ' DDCI T=-FF(TT 

EE ( I ) =0 • 

■“ "2" CONTINUE' ‘ • - ~ 

AA(N2)=HSUM/DX(N2)/16. 

CCT N 2 ) = RZ3 / 2 . / DX CN 2 1 F 

BB(N2)=-AA{N2)-CC(N2 ) 

— - FFTN2 } = FF('N2 1 ) ' ' 

DDCN2)=-ALAM*<H2+H1 ) /4. 

EE CN2T=DD~rN2~) +FFTN2 ) * - - — 

DO 3 I=NP1, NL1 

AACI )=T)13/2'*/DX ( IT ~ 

CC(I )=Hl3/2,/DX{ I— 1 } 

— 'bbt i fi-Asri r-ccn r — ' ' * “ ' 

FF ( I ) =Ai_AM*H 1/2 » 

D'cr{n=-FF{n 

EE( I) =0. 

3 CONTINUE 

AA(NP)=H13/2./DX(NP) ' 

CCTNP )=HSUM/DX ( N'2 ) /16~. ' * ~ 

BB ( NP ) =-AA ( NP ) -CC ( NP ) 

DD~( NP ) =TTDTNPT ) “ 

FF ( NP ) = ALAM* { H1+H2 ) / 4* 

” EECNP ) =DD ( NP ) +FF T NP") 

104 CONTINUE 

"" DO" 5 1=2, "Nil ----- - - — - 

J=I-1 

— X( -j } s'Axrn^p'oi 1+17 + dot n — ' ~ 

B l J ) -BB ( I ) *PO { I } +EE ( I } 

CTJT=CCTTr*P"0 (T-TT+F FIT) 

FI ( J)=A(J)*PO( 1+1 )+B( J >*POC I )+C( J)*PO( I— 1 ) 

— rnjT=~F'irJT 

5 CONTINUE 

— “t3'0”’6"TS-3TNL2T ' ~~ ~ “ 

J*I-1 

f roPT utxt =~ccrmE2'pFPu ( i-d+tfci ) • ■ 

F I DP ( J » 2 ) =BB ( I )*2.*PO( I ) + EE( I ) 

H — FTDPTJT3) = AAT IT*2T*P0TT+1 r+DDT I ) " “■ — 

6 CONTINUE 
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FIDPI l,2)=BB(2)*2.*PO( 2)+Et ( 2) 

FI DP ( 1 *3 >=AA ( 2 )*2.*PO( 3 )+DU ( 2 ) 

F I DP ( NL2 » 1 ) = CC (NLl )*2 **Pu( NLz >+FF (NLi ) 
FIDP(NL2»2) = bb(NLl)*2.*P0(NLl )+tElNLl) 
CALL TDLEQ ( F I DP , F I , DP , Nl2 ) 

DO 7 1=2, NLl 

PO( I )=PO( I )+DP( I — 1 ) 

IF (PO(I)-O.) lul, 101, 7 

101 PRINT 17 

DO 23 IP=1, NLl 

PRINT 18, IP,DX( IP) ,P0( IP) , up ( IP — 1 ) 

23 CONTINUE 
STOP 

7 CONTINUE 
RMAX = T0t-/2. 

DO 8 K = 1 » N L 2 
R=DP ( K ) 

IF (R-RMAX) 8, 8, 102 

102 RMAX=R 

8 CONTINUE 

IF (RMAX-TOL) 103, 103, 104 

103 CONTINUE 
PRINT 15 

DO 9 I P = 1 »NLl 

PRINT 16, IP, DX ( IP) , PO (IP) 

9 CONTINUE 
PRINT 19 
RETURN 
END 


SUBROUTINE TDlEO (A*C*X»N! 

SUBROUTINE SOLVcS A TRIDIAUUNAL bYSTti'i OF LINEAR LuUATIONS 
DIMENSION 4(50,3) , C < o 0 ) , X(5u), b(50,2), DO0) 

START 
J = N 

B l 1 * 1 ) = A ( 1 , 2 ) 

8(1,2) =A (1,3) 

D( 1 )=C( 1 ) 

JJ=J-1 
DO 5 K=2 , J J 

IF (ABS(B(K-1»1) )-ABS(B(K.-l,2) ) )3, h , 4 

3 B(K»l)=A(K*2)*B(K-l,l)/B(i^-l»2)-A(fc,l) 

B(K,2)=A(R,3)-tB(K.— 1,1)/8(K— 1,2) 

D( K) = (C C K J *b(u-l , 1 )-A ( K, 1 ) *0 (R-i ) )/o(!C-l ,2 ) 

60 TO 5 

4 B(K»l)=A(k»2)-A(K>l)*B(K-l*2)/B(N-l*l) 

B(K,2)=A(K,3) 

D(K)=CIK)-A(K,l)*Dlx-l)/bD.-l,i) 

5 CONTINUE 

7 X(J)=(C(J)*B(J-l,l)-A(J,l)*0(J-l))/(A(J,2)*b(J-l»l)-A(u,l)*blJ-l,2 
1 ) ) 



K=J-1 

15 X ( K ) - ( D ( K ) ~3 ( K » 2 ) *X ( K+ 1 ) ) / B ( K. » 1 ) 
IF l K-l I 100, 100, 10 
10 K=K-1 
GO TO 15 
100 RETURN 
END 
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Card 1 
31 

ALAM 

HI 

Tol 

NL 

N2 

Card 2 

VI(I) 

Card 3 

DX(I) 

Card 4 

P0(I) 


Format (5F10.4, 315) 

Ratio of B^/B. 

Bearing number, 

Normalized gas film thickness, 

Smallness for convergence. 

Total grid number for the step pad including two 
end points. 

Total grid number for the left edge of the step pad 
with H = H^. 

Format (8F10.5) 

Array of the squeeze numbers to be solved. 

Format (8F10.5) 

Array of increments of X defined as = ^i+1 "^i 

N _ 

such that 2 DX(I) = 1. 

1=1 

Format (8F10.5) 

Array of the guessed zero order pressure profile Po, 
on the pad with Po, 1 = Po, n = 1. 



